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Abstract 

Fluid flow applications can involve a number of coupled problems. One is the simula¬ 
tion of free-surface flows, which require the solution of a free-boundary problem. Within 
this problem, the governing equations of fluid flow are coupled with a domain deforma¬ 
tion approach. This work reviews five of those approaches: interface tracking using a 
boundary-conforming mesh and, in the interface capturing context, the level-set method, the 
volume-of-fluid method, particle methods, as well as the phase-field method. The history 
of each method is presented in combination with the most recent developments in the held. 
Particularly, the topics of extended finite elements (XFEM) and NURBS-based methods, 
such as Isogeometric Analysis (IGA), are addressed. For illustration purposes, two applica¬ 
tions have been chosen: two-phase flow involving drops or bubbles and sloshing tanks. The 
challenges of these applications, such as the geometrically correct representation of the free 
surface or the incorporation of surface tension forces, are discussed. 

Keywords: free-surface flow, interface capturing, interface tracking, NURBS, XFEM 
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1 Introduction 

In the numerical analysis of fluid flow, we often encounter free-boundary value problems: apart 
from the flow solution in the bulk domain, the position of (a portion of) the boundary is also 
unknown. This boundary can either be an external boundary or an interface between subdomains. 
At the boundary/interface, certain boundary conditions need to be fulfilled, which specify the 
position of the boundary. These conditions relate the variables of the flow (velocity, pressure, 
stress) across the domains under consideration of external influences, such as for example surface 
tension. 
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Figure 1: Google Ngram Viewer: quantitative analysis of popularity of five main interface de¬ 
scription methods. 


Numerically, in order to be able to compute the flow solution as well as the boundary/interface 
geometry, a measure to track the boundary starting from an initial position needs to be incorpo¬ 
rated. Recently, Coutinho gave an overview of the most common mesh-based methods in this 
area - level-set, volume-of-fluid, and phase-field - from a very different point of view: a pure 
quantitative analysis of popularity using Google books Ngram Viewer lISTI . Ngram Viewer is a 
tool, which analyzes all 30 million books digitalized by Google with respect to certain keywords 
II157II . The return value is the fraction of these books, in which the specific keyword appears. 
For this paper, we have added particle methods and interface tracking to the search, resulting 
in Figure [T] From the diagram, we get the notion that some methods already have a very long 
history of unremitting usage, such as the phase-field method, which - as we will later see - dates 
back already to 1871. The use of the level-set method rocketed within the last two decades, no 
doubt due to the large number of very effective recent developments. From pure numbers, it 
seems as if the volume-of-fluid methods may have already passed their zenith. Particle methods 
and interface tracking have so far not been able to compete with the more established methods on 
the grounds of popularity, but seem to be reserved for niche applications. Notwithstanding these 
crude quantifications, this paper aims to highlight the advantages and disadvantages of the above 
methods, illuminating the history of each method, but focusing on the most recent advancements 
in the individual fields. 

Numerous fluid flow applications involve deformable domains: drops and bubbles, die swell, 
dam break, liquid storage tanks, dendritic growth, spinodal decomposition, up to and including 
the topic of topology optimization, which can rely on the same boundary tracking methods. 

The scenario we consider here are flow solutions with either one or with two immiscible fluids. 
The general case is illustrated in Figure]^ Consider a cl-dimensional computational domain 
G C with boundary F = dO.. This domain contains one, two or more immiscible fluids, which 
are enclosed in subdomains G,(t). Position, number and shape of the individual subdomains may 
vary over time, i.e., both the domains themselves and the flow field are part of the solution. The 
main task in this context is to account for the interface P"', which separates the distinct fluid 
domains and is generally in motion. Apart from its exact position, the computation of geometrical 
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Figure 2: Illustration of the general two-phase flow scenario. 


quantities of as for example its normal and tangential vectors n and t or its curvature K are 


of interest. The corresponding computational tasks can be summarized as: (1) define the shape 
and location of the interface, (2) track the time advancement of the interface, and (3) set boundary 
conditions along the interface. 

The paper is organized as follows. In Section]^ the governing equations for the fluid flow 
with the appropriate boundary conditions for the interface are introduced. Particular challenges 
that come with the numerical treatment of these equations are highlighted in Section The 
subsequent sections concentrate on the five numerical methods considered within the scope of 
this paper: particle methods in Section]^ the volume-of-fluid method in Section|^ the level-set 
method in Section|^ the phase-field method in Section]^ and mesh-conforming interface tracking 
in Section]^ To illustrate the methods, two common applications, drops and sloshing tanks, are 
the topics of Sections and [T0| 

2 Governing equations 

In general, the governing equations for incompressible fluid flow are the instationary, incom¬ 
pressible Navier-Stokes equations. Consider an instationary fluid flow problem with any number 
(in this paper, one or two) of immiscible Newtonian phases. The computational domain at each 
instant in time, denoted by , is a subset of where nsd is the number of space dimensions. 
Then at each point in time f G [0, T], the velocity, u(x,t), and the pressure, p(x,f), in each phase 
are governed by the following equations: 



( 1 ) 


V-u = 0 on (fl,); Vf G [0,r], 


( 2 ) 


for i = I,... ,np (number of phases) and p,- as the density of the respective fluid. The stress tensor 
Gi is defined as 
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( 3 ) 


ai{u,p) =-pl + 2Hie{u) on (Hf),-, 

with 

e(u) = l(Vu + (Vu)^), (4) 

where /r, denotes the dynamic viscosity, f includes all external body forces referred to the unit 
mass of fluid. The computational domain is divided into parts (flf),, each occupied by fluid /. 
Note that the spatial domain is time-dependent, which is indicated by subscript t. A subdomain 
may never be occupied by more than one fluid. The interface of two subdomains, dQ.i n dQ. 2 , is 
denoted by r“. 

To allow a better grasp of the Navier-Stokes equations we will illustrate the purpose 

of the individual terms, following in large parts the description in |[46l . 

• Temporal derivative ^ : The change in velocity with respect to time. This change over 
time is governed by the following influence factors: 

• Inertia term u-Vu : This term is a convection term arising from the conservation of 
momentum. The momentum of each portion of fluid needs to be conserved. Therefore, it 
needs to move with the fluid - it is convected with the fluid. 

• Pressure term —Vp : This term appears when the Newtonian constitutive equation © is 
inserted into Equation ([T]). It includes forces resulting from pressure differences within 
the fluid into the formulation. Especially in the case of incompressible fluids, it is very 
important to consolidate this term with Equation ©. 

• Friction term /r,V^u : Again, this term arises when the Newtonian constitutive equation 
© is inserted into Equation ©. We obtain a diffusion operator equalizing the velocity of 
neighbouring elements. The more viscous the fluid, the stronger is the friction between 
neighbouring particles and thus the equalizing effect. 

• External forces f : This term includes external body forces such as gravity or vibrational 
excitation. 

Eor creeping flows (i.e., Reynolds number <C 1), the advective term in the Navier-Stokes 
equations is often neglected, giving rise to the Stokes equations. If furthermore the solution 
remains unaltered over time, the stationary Stokes equations can be employed: 


(5) 

( 6 ) 

The constitutive equations for Newtonian fluids remain the same as above. 


—V<T,=f on D.i, 
V • u = 0 on fl,. 
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2.1 Boundary, initial, and interface conditions 

In the transient case, a divergence-free velocity field for the whole computational domain is 
needed as an initial condition: 


u(x,0)=u°(x) inflatt = 0. (7) 

In order to obtain a well-posed system, boundary conditions have to be imposed on the external 
boundary of Q., denoted as F. Here, we distinguish between Dirichlet and Neumann boundary 
conditions given by: 


u = u onr„, t e [o,r], ( 8 ) 

n-<T = h on F/j, t G [0,r], (9) 

where u and h are prescribed velocity and stress values. F„ and F/, denote the Dirichlet and 
Neumann part of the boundary, forming a complementary subset of F, i.e., F„ UF/, = F and 
Fm n F/, = 0. n refers to the outer normal vector on F/,. 

In the case of two-phase flow, at the interface between the two phases, we impose interface 
conditions, which couple the velocity and stress between the two domains. The first interface 
condition implies that the velocities are continuous across the interface: 


Ui =U2. 


( 10 ) 


The second interface condition is based on the Laplace-Young equation |[T^ : 


(ff2-ffi)ni = yx-ni. (11) 

Here, ni is the, with respect to Di, outward unit normal vector on FJ”^ 7 the surface tension 
coefficient and K the sum of the principal curvatures of Ff"h 

In addition to the interface conditions, the jump in density and viscosity between the fields 
leads to non-smooth behavior of the field quantities. If we examine condition ( fTO] ) closer, we 
note that it states that the normal component of the velocity = u • ni as well as the tangential 
velocity component Uf = u- ti must be continuous across the interface. Unaffected by this, it can 
however be shown that the normal gradient ^ of the tangential velocity is discontinuous II137II : 


dufi duf2 
dn dn 


(Ml 



( 12 ) 


with denoting the tangential derivative. In consequence, the velocity has a kink across the 
interface if /ti / /i 2 - 

Inserting the definitions for stress and strain (Equations Q and Q) into interface condition 
( [TT] ) results in: 
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[-pil + Pi (Vui + (Vui)^) + P2I - At 2 (Vu 2 + (Vu 2 )^)]ni = yKiii. 


(13) 


If we restrict ourselves to only the normal component of relation ( [131 ), an expression for the 
pressure jump across the interface can he derived: 


n[[-piI + /ii(Vui + (Vui)^)+P 2 l-At 2 (Vu 2 + (Vu 2 )^)] -ni = yfc, (14) 

44>-pi+P2 + [2iUi(Vuini)-ni-2^2(Vu2ni)-ni] = YK, (15) 

p d Ufi 1 d Ufi 21 

<^-pi+P2 + [2pi-^-2p2-^] = yK. (16) 

From Equation ( [T^ we can deduce that the pressure jump across the interface depends on the 
surface tension coefficient, the curvature of the interface, and the jump in the viscosity weighted 
hy the normal derivative of the normal velocity component |[T^I117[I1241I . In particular, this 
has the consequence that pressure jumps across the interface are possible even when no surface 
tension effects are present, a possibility often disregarded in the modelling process. Usually 
however, this effect is indeed negligible as either the difference in viscosity or the change in the 
normal velocity are not high enough to actually make an impact II117II . 

A further influence on the pressure distribution is the volume force per unit mass of fluid f. To 
see this, we examine the hydrostatic case, i.e., u = 0, where the momentum equation o reduces 
to 


Vp; = p;f, infl. 


(17) 


Focusing on the interface a jump condition for the pressure gradient is obtained: 


Vpi — Vp2 — Pifi — P2f2 — (Pi — P2)f- (18) 

The last simplification can be made since the volume force is usually constant across the entire 
domain. We note that we obtain a jump in the pressure gradient proportional to the jump in 
density, implying a kink in the pressure distribution. This behavior is extendible to u / 0. 

All effects combined (depicted in Figurej^, in the general two-phase flow scenario we have to 
account for a kink in the velocity as well as a jump and kink in the pressure at the interface. 

2.2 Variational form 

In the finite element method, we work with the variational form of the governing equations. 
One important measure when deriving the variational form is the integration by parts with 
the goal of reducing the differentiability requirements on the trial functions (i.e., the solution). 
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Figure 3: General situation in two-phase flow: A kink in the velocity (due to fXi ^ H 2 ) and a jump 
(due to surface tension) as well as a kink (due to volume forces) in the pressure across 
the interface. 

Integration by parts furthermore naturally includes Neumann type boundary conditions as in our 
case Equations @ and ( [TT] ). II113II 

The variational form of Equations ([T]l-(|^ can be expressed as follows: Eind u and p such that 


Vw, yq: 



(19) 


Sections |^-[^ will introduce a variety of methods that can be employed to solve the above 
equation numerically. 

3 The challenges in free-surface flow 

In this section, we will analyze particular challenges of Equation ( [T9l ): the free boundary problem 
it represents, the treatment of the surface tension term, as well as the solution of additional 
equations on the interface. 

3.1 The free boundary problem 

Equations ([T]l - Q describe a scenario where a partial differential equation is solved for an 
unknown function - in this case the fluid velocity u(x,t) and the pressure p{x,t) - but at the 
same time the exact extent of the computational domain (or portions of it) is also unknown. 
We are dealing with a free boundary problem. Consequently, one of the computational tasks 
will be to monitor the domain shape - and especially the interface/boundary - throughout the 
simulation. Broadly, two significant approaches can be distinguished: interface capturing and 
interface tracking. The main difference can be located in the viewpoint, which can be either 
Eulerian or Eagrangian, as we will see in the following two sections. 
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3.1.1 Interface Capturing 

Interface capturing approaches are based on an Eulerian description and widely used in deformable 
domain problems. They define the interface implicitly on a fixed mesh. A characteristic scalar 
field <p is used to identify the two phases as well as the interface along the boundaries of the 
individual fluid domains. Depending on the different methods, this scalar field may be described 
for example by a discontinuous Heaviside function or a signed-distance function. In order to 
account for the interface motion, a standard advection equation 

d(j) 

—+u•¥(/)= 0 inll, te 0,r , (20) 

at 

is solved with u as the fluid velocity. The most common representatives of this category are 
particle methods (Section]^, the volume-of-fluid method (Sectionj^, and the level-set method 
(Section]^. In addition, the phase-field method (Section[7]l shares some features with the interface 
capturing methods. 

A great advantage of the interface capturing approaches is that they are inherently able to 
account for topological changes of the interface. This allows for a much more flexible interface 
description than in interface tracking approaches. What needs to be considered however are the 
aspects of discontinuity treatment across the interface, mass conservation, and the application of 
boundary conditions along the interface, which remain a challenge in interface capturing. 

3.1.2 Interface Tracking 

In interface tracking approaches, the interface is described explicitly on a boundary/interface 
conforming mesh. The idea is to track the position of the mesh nodes x in a Lagrangian fashion 
by integrating the evolution equation 


-^=u(x,t), (21) 

with u as the fluid velocity in the domain. We can distinguish between fully Lagrangian ap¬ 
proaches, where all mesh nodes are treated in a Lagrangian fashion, and Arbitrary Lagrangian- 
Eulerian (ALE) approaches, where the Lagrangian treatment is only applied to a portion of the 
mesh nodes, usually those along the boundary/interface (cf. Sectionj^. 

Interface tracking approaches offer an accurate and computationally efficient approximation of 
the boundary/interface. Furthermore, the imposition of boundary conditions at the interface is 
simple with mesh nodes lying on the interface itself. However, firstly the mesh quality will usually 
degrade significantly in the course of large deformations and secondly topological changes require 
special treatment. In both cases, the last resort is remeshing with respect to the new interface. 
One consequence of remeshing is the need to project the field values from the old to the new 
mesh, a ceaseless cause for errors, and at the same time computationally costly especially in 3D. 
Consequently, it is desirable to keep the frequency of remeshing low (even to the level of no 
remeshing) II208I . 
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3.2 Surface tension 

Much of the computational accuracy and stability of two-phase flow depends on the discretization 
of the term responsible for the surface tension force; in the variational form Equation ( [T^ this is: 


7 / KVf-ndx, ( 22 ) 

Jr;'” 

with surface tension coefficient 7 , K as the sum over all principle curvatures, w denoting the test 
function and n the normal vector. For most numerical approaches - although we will see some 
exceptions later - the evaluation of the curvature K is nonviable as it contains second derivatives. 
As an alternative, two main strategies have evolved to avoid the direct evaluation of Equation ( |22| ): 
the Laplace-Beltrami technique and the Continuum Surface Force approach. These approaches 
are applicable irrespective of the interface discretization method. They will be detailed in the 
following. 

3.2.1 The Laplace-Beltrami technique 

The idea of the Laplace-Beltrami technique goes back to Dziuk ll 6 ^ and has been closely 
described in lIT^ ISTl 1^ . The method replaces the above surface integral ( |22| ) requiring the 
curvature (and therefore the existence of second derivatives) by an integral requiring only gra¬ 
dient computations. In the following description, we will show two derivations of the method. 
The first is relevant when after discretization, an interface-conforming mesh (explicit interface 
representation) will be used and the second is suitable for implicit interface descriptions. 

When the interface will eventually be resolved explicitly by the mesh, our starting point is the 
following relation regarding the curvature vector K of the interface |[73l : 


k = A^X, (23) 

X:DCM“^^^ (24) 

where X is an immersion, i.e., a mapping from an arbitrary parameter space to a boundary 
or interface. An immersion can be understood as a differentiable map between differentiable 
manifolds whose derivative has full rank. In the case of free-surface flows, we deal with smooth 
embeddings, which are injective immersions. The curvature vector K always points in normal 
direction: It can therefore be written as (ten). We obtain: 


fcn = A«X. (25) 

This means that we have arrived at the expression needed in term ( [22l ). Again in the case 
of smooth embeddings, A^ takes on the form of the so-called Laplace-Beltrami operator, a 
generalization of the Laplace operator to manifolds, whose exact definition will be given later in 
Equation ( |2^ . Before doing so, we define a quantity essential for any geometrical computation 
on our immersion X: the metric g given as 
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(26) 


g = (gij) = JX^JX e 

The metric is intimately related to the measurement of distances and angles, and as such, also 
an important requisite when integrals are to be defined in such a way that they are invariant to 
coordinate transformations. For an arbitrary scalar function /, the coordinate transformation 
yields: 


/ f{x)dx‘^= / f{X{^))\l Ae.ig{^)d^. 
IV"' JTf 


(27) 


Note that in the cases, where the coordinate transformation maps W‘^‘‘ i—)■ ]R”“, JX is square. 
As a consequence, det(g) = (det(7X))^ and the definition includes the well-known transformation 
(change of variables) rule used routinely within the finite element method. 


Within the Laplace-Beltrami technique, the general procedure in modifying the integral in ( |22| ) 
will now be to replace fcn by the Laplace-Beltrami operator applied to the immersion X, integrate 
by parts and then transform the integral to reference coordinates. The definition used here for 
the Laplace-Beltrami operator can be found in fT^ . The Laplace-Beltrami operator applied to a 
scalar function / is: 


djf), 

with (g’j) = Inserted into the surface integral and integrated by parts we obtain: 


(28) 



ffw -n dx 


(^) 


7 [ w-A^Xdx 

Jpm 


j w di {g‘^ y/d^ djX) dl 
JV' 


B 

B 





(29) 


The steps in Equation (22 1 are detailed in the following: 


(1) The curvature vector is replaced using the Laplace-Beltrami teehnique. 

(2) We insert definition ( [2^ for the Laplace-Beltrami operator. The int^ral is transformed to a 
parametrization of the interface in ^-coordinates as given in Figure |4| This parametrization 
is still smooth. 


(3) Under the assumption that w is sufficiently smooth, integration by parts can be performed. 
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Figure 4: The interface P"', which is given in x-coordinates, is reparametrized with the parameter 
^ G [0,1]. This is a smooth embedding and the Laplace-Beltrami technique can be 
applied. Partial integration can be performed. The integrals are then as usual computed 
on the reference element in ^ coordinates. The ID reference element is mapped to a 
portion of the boundary of a 2D mesh (here, we see an excerpt). 


(4) The interface is discretized using a boundary conforming finite element mesh. The integral 
is computed element-wise through a reparametrization into element reference coordinates 

For the integration by parts, we assume that the surface is closed and therefore, the boundary 
term vanishes. This is an assumption regularly made lISTlIWl . 

In order to give an example of the actual implementation, we will demonstrate how to evaluate 
the immersion term in Equation ( |29l ) in a standard hnite element setting. is defined as the 
mapping from reference coordinates ^ to physical coordinates x: 


X^:^^x X'^(^) = £iV,(^)x^. (30) 

k=i 

Here, denotes the finite-element shape function for node k and nen indicates the number of 
nodes per element. Note that we refer to the boundary element as indicated in Figure 

The derivative of X with respect to the (here only one) reference coordinate ^ evaluated at a 
specific point is then: 


5X 


dNk 

h 


Xk. 


(31) 


This concludes the explicit case. 

In the implicit case, the starting point is to express the sum of principle curvatures of the 
interface F'”' by a relation similar to Equation ( |25] ) as Wi\ IFtI : 
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Figure 5: The tangential gradient of a function / can be evaluated by computing the full gradient 
and subsequently subtracting the normal component of the gradient. 


ten — . (32) 

In Equation ( [^ , idpm signifies the identity mapping on T"'* defined as idpm (x) : Y'"’ i—)• F'”^ := 
X = (xi,X 2 , ■ ■ ■ ,Xnsd)^. Note that this definition implies that V {idpm) are exactly the basis vectors 
of Furthermore, Ap denotes again the Laplace-Beltrami operator. However in this variant, a 
different definition is employed. This is also the reason, why for now a different symbol has been 
chosen. Apart from the definition in Equation ( [2^ , it can also be defined using the tangential 
derivative and the tangential divergence ll97ll : 


Ap = Vp-Vp, (33) 

where the tangential derivative Vp of an arbitrary function / is defined as: 

Vr/:=V/-(n-V/)n=(I-nn^)V/, (34) 

and the tangential divergence of f is defined as : 


Vp-f := V-f-n^(Vf)n. 


(35) 


Using these relations to replace ten, we can again take advantage of the idea that integration by 
parts within the weak form can be used in order to reduce the order of derivatives. Again, the 


integral in equation ( [221 ) is restricted to the interface and therefore the integration by parts is not 
straightforward. In the case of a closed surface where the boundary term vanishes, this leads 
to 


7 / Kw-ndx = —Y Vpw : Vp/r/pm, r/x = — 7 / tr((I —nn^)Vpw) r/x. (36) 

Jpm Jpm Jpm 
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The assumption of a closed surface is valid in many applications where surface tension 
is important, such as for example drops and bubbles. Consult ll97l for derivations for open 
boundaries and modifications for moving interfaces. Note that as compared to the derivation 
in the explicit case, the extraction of the tangential gradient using (I — nn^) is necessary. This 
stands contrary to the definition in Equation ( [29l ), where all gradients already point in tangential 
direction. 

The obvious question to pose at this point is how the two presented expressions for the Laplace- 

1), can be related. In the context of 


Beltrami operator, in Equation pS] ) and Ap in Equation i 
this method, the crucial difference between the explicit and the implicit approach is that in the 
explicit case, the test function w can be expressed directly in terms of the reference coordinates 
^; even along the interface. In the implicit case, however, the test function on the interface w(^) 
needs to be expressed in terms of the test function in the bulk domain, which we will denote as 
w(x) : 


w(^) =w(X(^)), 
w : X(D) C : 


(37) 

(38) 


Over the course of the next equations, we will demonstrate that despite this difference be¬ 


tween the approaches. Equations ( [29| ) and ( [3^ are to be considered as equivalent. If we insert 
definition (|T7|) into Equation ([29l), which has so far only been used for the explicit case, we obtain 




dx^ 

-w* ) o‘j 


d ...Y 5 d 




(P(«))) 

= (vw|xq^)P(^)) y/detg 

Based on both the definition of the normal n and {g‘^) = 


(39) 


and 


Pn = 0, 


p ^ ye _ ^ -ye 


(40) 


(41) 


hold. In other words, P is a projection onto the tangent space. Since in the 2D case {n, ^X'^} 


and in 3D {n, g^X'^, ^X*^} form a basis. Equations (^l - ( |4T] ) define P unambiguously. In 
particular, this means that 
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P = I - nn^ . 


(42) 


With Equation ( |3^ and = P we obtain 

—7 [ tr(PVrw) \/detg = —7 / tr(PVrw)(ir. 
j/ Jr'" 

The latter is equivalent to Equation ( |3^ . 

Note that in all cases of both explicit and implicit boundary descriptions, the assumption that 
Equations ( [22l ) and either ( [29l ) or ( |3^ are completely equivalent only holds if is sufficiently 
smooth. This will in general not hold for the discretization of E"" WT\ . 

3.2.2 Continuum Surface Force approach 

Another alternative to Equation ( [2^ is the Continuum Surface Force approach (CSF) as introduced 
in |[32l . This is of particular interest when modelling topologically complex interfaces. The basic 
idea is to replace the interface coupling condition ( [TT] ) by a localized volume force term in the 
momentum equation ([T]l. The force term has the following form Il96l : 

fr = ytcSrn. (44) 

Here, 5r refers to a smoothed Dirac 5-function needed to select the surface force to the 
proximity of the interface, i.e., it selects a narrow band around the interface, fr is designed in 
such a way that it exactly replicates the surface tension force per interfacial area that the surface 
integral would generate. Away from the interfacial region, the surface force is 0. Throughout the 
transition region, all fluid quantities vary continuously. 

3.3 Solving additional equations on the interface 

In some applications, it may be necessary to solve additional partial differential equations on the 
interface/boundary: an example for the topic of solving partial differential equations on hyper¬ 
surfaces (manifolds). One prominent example for such a case is the solution of the transport 
equation (advection-diffusion equation) as an additional equation to the governing equations 
of fluid mechanics lIMl II191 [14011 . materials science |l5ll|76j|, and biology ifTTl 11341 [14511 . The 
main topic is the consideration of insoluble surfactants (surface active agents), which adhere to 
the phase interface and influence the fluid flow in an either restraining or stimulating way. One 
type of surfactant known to all of us in our everyday lives is soap. It is important to note, that 
the surfactant is located exclusively at the interface/boundary and never in the bulk domain (cf. 
Figure]^. Surfactants appear either as pollutants, which have accidentally entered the flow, or are 
purposefully added for example to decrease the surface tension and thereby increase the ability of 
a liquid to wet solid surfaces 15^ . 
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Figure 6 : Depending on the flow field and the diffusion coefficient, the surfactant concentration 
on the interface of, e.g., a drop, can he determined in every time step. 


The two main driving forces for the transport are advection (movement of the interface) and 
diffusion (molecular diffusion along the interface) WT\ . Even though the governing equations for 
the scalar transport processes are often known and well understood, solving the respective partial 
differential equation on arbitrary and even moving manifolds is not straightforward. The effects 
brought about by the embedding, as already mentioned in Section [3. 2.1 [ have to be accounted 
for. The solution approaches differ significantly depending on whether an explicit or implicit 
mesh description has been chosen (cf. Section [3T| ). In the explicit approach, an interface mesh 
can be extracted from the bulk mesh thus easily defining the domain on which the transport 
equation is to be solved. In order to solve the proper equation, the differential operators of the 
partial differential equation have to be generalized in order to account for the arbitrarily shaped 
hyper-surface ll64llMl . Note that during the integration on the reference domain, it is mandatory 
to include the metric introduced in Equation ( [27] ). Implicit approaches ||2| [^ [65] 11691 122 111 
typically extend the scalar field away from the interface to the whole domain in order for all 
differential operators to be defined properly. One possible approach to the extension can be found 
in ll3^l4T1l . This obviously increases the dimension of the problem, but the transport equation 
can then be solved in the whole fixed computational domain. llWl provides a good overview of 
different Eulerian and Lagrangian approaches. 

The solution of the transport equation requires a coupling of the flow solution with the transport. 
In addition, in many applications we would like to consider the effects that the surfactant has 
on the interface; i.e., a coupling of the transport solution to the flow equations. This coupling 
enters through the surface tension coefficient 7 , which can be considered to be a function of 
the surfactant concentration G. This leads to tangential forces along the interface, the so-called 
Marangoni convection. Il83l l97ll 

The standard advection-diffusion equation would be given as: 

^+V-(aG)-V-(dVG) = 0, (45) 

where a denotes the (possibly time dependent) advection velocity and d the diffusion coefficient 
matrix. In the case, where the above equation needs to be solved on a deformable hypersurface, 
some modifications need to be taken into account ||83l |97ll : 
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1. In the surfactant context, we will usually assume that d is constant and isotropic. Therefore 
V • (dVG) = dAG. 

2. Since the equation is to be solved on an interface embedded into a larger domain, the 
differential operators in Equation ( |45] ) need to be replaced by the corresponding operators 
restricted to the interface: the gradient V is replaced by the tangential gradient Vp (cf. 
Equation (|3^), the divergence V- by the tangential divergence Vp- (cf. Equation ( |35| )) and 
the Laplace operator A by the Laplace-Beltrami operator Ap (cf. Equations (^l and ([^). 


3. The advection velocity - and in the case of surfactant transport the advection velocity will 
be the fluid velocity u - needs to be restricted to its tangential component up = (I — nn^)u. 


4. The term Vp • (uG) can be rewritten as 


Vp • (uG) = GVp • u + Up • VpG . (46) 

In the cases of deformable interfaces, the first term ensures conservation of mass even 
under local changes in the free-surface area. Note that it is important to use u instead of up. 
Otherwise, a change in surface area would not be possible. If the interface is fixed, this 
term is 0 as u • n = 0. 

This leads to the final equation for the surfactant concentration G: 

+ (up • Vp)G — ApG + GVp -0 = 0. (47) 

The above equation may be rewritten in many ways l|2j|64l|83l, making it more accessible to 
the different numerical schemes. 


4 Particle Methods 

Particle methods utilize mass-less particles distributed in an Eulerian mesh to capture the fluid 
flow and in particular the interface position. The most prominent particle method is the Marker- 
And-Cell (MAC) method developed by Harlow and Welch II102I . It is one of the first examples of 
an interface capturing method; the numerous examples presented in |102|| attest to the method’s 
extraordinary flexibility. 

The original approach was developed for a single phase - i.e., / = 1 in Equations ([T]l - Q - 
thus modelling a fluid surrounded by air. In later versions, the method could be extended to multi¬ 
phase flow. In all cases, the particle methods solve the “one-fluid” version of the Navier-Stokes 
equations, where the surface tension force is included directly and the fluid properties, e.g., p and 
p, vary according to an indicator function. 

The method covers with a computational mesh the maximal possible fluid domain and subse¬ 
quently distributes marker particles throughout the entire domain currently covered with fluid (cf. 
Eigure|7]). Any cell not containing marker particles is identified as an empty cell. Cells with at 
least one marker particle and at least one common boundary with an empty cell are the interface 
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Figure 7: In the Marker-and-Cell method by Harlow and Welch II102L the following convention 
is used for the identification of the fluid domain and the empty domain: Any cell not 
containing marker particles (blue dots) is identified as an empty cell. Cells with at 
least one marker particle and at least one common boundary with an empty cell are the 
interface cells (marked in red). Cells accommodating at least one marker particle and 
furthermore surrounded only by other cells containing marker particles are marked as 
fluid cells. 


cells. As a last category, cells accommodating at least one marker particle and furthermore 
surrounded only by other cells containing marker particles are marked as fluid cells. 

Throughout the simulation, the marker particles are advanced with the local fluid velocity 
interpolated from the Eulerian grid. 

The MAC method is strongly connected to finite-difference discretizations. In principle, 
however, it can also be combined with other methods, such as the finite element method in 1901. 

In the MAC method, the velocity and pressure solutions are usually obtained sequentially. First, 
a preliminary velocity field that we will denote with 0 "*“” is computed including all effects of 
the momentum equation (Equation Q) except for the pressure: Of the terms listed in Section]^ 
the inertia term, the friction term and the external forces contribute. The relevant equation reads 
(adapted from [T] and Bfiil '): 


— = -(u”-V>"+f+ vAV. (48) 

Here, the superscript h on the operators indicate that they are discretized in an appropriate 
fashion, e.g., with finite-difference schemes. The superscript n denotes the value of the corre¬ 
sponding field at the previous time step. Note that Equation ( [4^ can either be solved in one step 
or subdivided into several steps, computing the contributions of the individual terms one by one 
as done in ||46l. 

The computed u'"'’'” will most likely not be divergence free and therefore. Equation Q is not 
fulfilled yet. However, the contribution of the pressure term to the velocity field has not been 
considered so far. This contribution is still available to adjust the velocity field by setting the 
pressure values appropriately. and its correction factor — Vp are inserted into Equation Q: 
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(49) 


Again, the superscript h refers to a discretized operator. Equation ( [491 ) is a Poisson-type 
equation for the pressure. By solving this equation, the pressure values and with those, 
can he obtained. This procedure is comparahle to projection methods - or Chorin-Temam method 
- utilized in the context of the finite element method ( l[74l and references therein). 

Notable is the development of SMAC (simplified-MAC) 161 only five years affer fhe original 
MAC implemenfafion, which reduces fhe complexify of fhe original approach parficularly when 
if comes fo incorporafing differenf fypes of boundary conditions. Ifs secref lies in fhe compufafion 
of fhe pressure field nof based on a preliminary velocify field, buf rafher a pofenfial funcfion. 


4.1 Properties 

Wifh regard fo fhe fypical applications in fluid flow, MAC has fhe advanfage fhaf fhe fracer 
particles are indisfinguishable from one anofher: The joining and separation of fluid parfs can 
be effortlessly simulafed. Neverfheless, if suffers from drawbacks especially regarding fhe 
descripfion of free surfaces. Sfress effecfs af fhe free surface, such as surface fension and confacf 
angles are difficulf fo include. This is parficularly frue since, if used wifh finife difference, MAC 
is enfwined wifh a discretization on sfaggered grids, i.e., storing velocify on cell faces and fhe 
pressure af fhe cell cenfer, which has been found fo provide higher sfabilify ll46l . Boundary 
conditions af fhe free-surface can fhen for example only be sef as pressures af fhe cell cenfer of 
all cells identified as inferface cells ll66ll . i.e., nof af fhe acfual boundary, or fhey may require 
complex surface tiffing fechniques such as fhose infroduced in li55l . In addifion, if is difficulf fo 
compufe normals and curvafure from fhe parficles II138II and fo confrol particle driffing II224II . In 
general, fhe mefhod can suffer from sfabilify problems, which may be difficulf fo defecf, wifh 
apparenfly reasonable particle disfribufions resulting from gross approximation errors. Some 
modification were proposed fo sfabilize fhe MAC mefhod in li34ll . McKee HI541 identifies fwo 
furfher consequenfial resfricfions: MAC could nof be applied fo arbifrary domain shapes and was, 
for a long time, resfricfed fo fwo space dimensions. The laffer was due fo lacking efficiency in fhe 
solution mefhods for fhe correcfed pressure, preventing fhe use of a large number of parficles. In 
addifion, care had fo be exerfed on fhe choice of fime-sfepping (no marker is allowed fo cross 
more fhan one cell per fime-sfep) 111541 . 


4.2 Recent advances 

Recalling fhe drawbacks lisfed in fhe lasf secfion, fwo main poinfs have recenfly been addressed 
in fhe area of MAC: fhe reslricfion to simple shapes in 2D and fhe topic of applying boundary 
conditions af walls and on fhe free surface. 

A developmenf in fhe area of fhe former are GENSMAC 12091 and GENSMAC3D 12101 . 
This extension provides fhe applicabilify of MAC to domains of arbifrary shape and dimension 
as opposed fo fhe sfraighf lines required in fhe original implemenfafion. This approach has 
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been applied to viscoelastic extrudate swell in II162112111121211 . i.e., an application where the 
components of the stress tensor appear as additional unknowns in addition to the only primary 
variables of velocity and pressure. HI9011 proposes an extension to an arbitrary number of phases. 
In a reduction of the mesh size and thus computational effort is proposed. The procedure is 
illustrated in Figure 

Boundary conditions on interfaces remained a challenge to MAC for decades. The only way 
to overcome this challenge was so far the combination with other methods, such as Lagrangian 
meshes (cf. Sectionj^, level-set methods (cf. Sectionj^ or volume-of-fluid (cf. Sectionj^ II154II . 
Santos Hi901 succeeded by including an additional Lagrangian mesh to track the free surface. ifTOl 
works with an ordered list of connected surface markers, which are advected along the streamlines 
of the flow field using a Runge-Kutta integration. Since the flow field will most likely disturb 
the homogeneity of the marker distribution (concentration of streamlines or vortical flow), they 
are subsequently added, deleted and redistributed. H22411 proposes a hybrid approach in such a 
sense, that the standard MAC method with its usual Eulerian grid is used for the fluid domain. In 
addition, special marker particles identify the interface, which is subsequently reconstructed using 
a Lagrangian mesh. Surface tension was also considered in ni93H . where the authors compared 
an implementation using a Lagrangian mesh with an implementation based on a particle level-set 
method. In H138I . a level-set approach is combined with Lagrangian marker particles as a means 
of level-set reinitialization. ni35i 11361 present a method where the interface is tracked using 
mesh-less particles. The particles are explicitly connected to certain points of an underlying 
Eulerian reference grid. Connectivity among the particles is not required however. In each 
step, the interface is resampled based on a local reconstruction of the interface, on the one hand 
redistributing the particles and on the other hand updating the connecting points in the reference 
grid. The approach taken in H218II is based on the point-set method first introduced in H2131 . It 
uses an indicator held (a function which is either 0 or I depending on the corresponding fluid 
domain), which is constructed from massless marker particles. Subsequently, the normals and 
curvatures can be computed from the indicator function, which has been smoothed using the 
reproducing kernel particle method (RKPM). A completely different approach has been adopted 
in ni32ll . where all gradient terms in the Navier-Stokes equations are computed based on the 
interaction between all particles within a certain kernel. This method has been extended with 
more efficient solvers in H123II . 

The MAC-method served as inspiration for the volume-of-fluid method described in the next 
section. 


5 The Volume of Fluid Method 

The volume-of-fluid (VOF) method was developed in 1981 by Hirt Hl05l . At the time, he worked 
in an environment where there were two options available to represent surfaces: either marker 
particles (cf. Section or height functions combined with line segments. Height functions 
combined with line segments - Hirt gives H108I and HI641 as references - were the first attempts 
to represent the actual interface. If the height function is used, the height of the interface as 
compared to a given reference line is stored for each discrete point in the domain. Problems occur 
for multiple-valued height functions as they would be seen for drops and when the interface slope 
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(a) Maximal domain 


(b) Reduced domain 


Figure 8: (a) In the original MAC approach, the maximal possible fluid domain (black cells) was 
meshed, even if only a minor portion of cells (red cells) contain fluid, (b) In more recent 
approaches, the mesh consists only of the cells containing fluid plus an additional layer 
of ghost cells. This is sufficient as there is a time-step requirement, which guarantees 
that fluid particles cannot move through more than one cell per time step. l!4^ 



Figure 9: Excerpt of a domain containing nine cells: The interface is indicated in blue. The area 
left of the interface is filled with fluid. In the VOF method, for each cell, the fraction of 
fluid contained in the cell is stored. Cells that are completely filled with fluid have the 
value 1. Cells without any fluid have the value 0. 


exceeds the grid resolution. The height function approach can be generalized to ordered chains 
of line segments. Although the two given issues with the height function are resolved now, this 
approach has its drawbacks when it comes to curve intersections or an extension to 3D. 

As described in the previous section, the MAC method does not define the interface, but rather 
the fluid regions, thus leading to no logical problems during interface intersection as well as a 
straightforward extension to 3D - apart from the excessive storage requirements, which were a 
concern to Hirt. Building on this scenario, he proposed the VOF method. This successor of the 
MAC concept replaces the discrete marker particles with a global field defining the location of 
the fluid(s). This global field is called the volume-fraction field F, sharply dividing fluid (f = 1) 
and empty {F = 0) areas. For each cell, F takes on one specific value 0 < f < 1, indicating 
the amount of fluid contained in each cell (cf. Figure]^. One way of defining F is the use of 
a characteristic scalar field 0. In the case of the VOF method, (/) = with as the 
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heaviside function. The volume fraction for each cell Q. can then be defined as: 


F(a,i) = ^ / Jif{x,i)da. 


(50) 


With only one value per cell, the storage requirements for this method are significantly reduced 
as compared to the MAC approach. A very important feature of the VOF method is its inherent 
mass conservation |[3^ . 

The advancement of the interface in time is governed by the following advection-type equation: 



(51) 


Due to its discontinuous character, the volume fraction function cannot be advected directly. 
Instead, it is connected with the mass conversation equation of the Navier-Stokes equations, thus 
including the volume-fraction field info fhe flow solufion. The resulfing equation can be reduced 
fo fhe compufafion of fluxes across fhe elemenf faces. In ifs mosf basic version, fhe VOF mefhod 
has fhree major difficulfies: 

1 . fhe compulation of shape derivatives (as for example needed for curvalure evaluation) is 
obslrucfed by fhe use of a non-smoolh inlerface represenfalion, 

2 . as fhe volume fraction function is advecled, fhe inlerface becomes more and more diffuse, 

3. imposition of boundary conditions along fhe inlerface is impossible Ill05[ r7]l. 

One possible solution to all three problems is the reconstruction of a sharp interface in each 
time step. Options for this reconstruction will be discussed in the following section. Nevertheless, 
the VOF method is still utilized today in its original version, as for example in Il44l . 

5.1 Recent Methods for Interface Reconstruction 

Many advances in the VOF method are based on the aspect of reconstruction of the interface line. 
Over the past decades, numerous approaches have evolved. This large number of approaches 
already hints towards one major problem: As the interface shape is inferred from the volume data 
it is never unique. The reconstructed shapes are in general either piecewise linear or piecewise 
constant. il88ll 

In the original implementation of VOF, a method that would later become known as SLIC 
(Simple Linear Interface Calculation) II165I1 was proposed, where the interface is reconstructed 
using line segments, which can be either parallel or perpendicular to the major flow axis. The 
direction of the line segment is defined through the coordinate direction in which a larger change 
in the fluid volume occurs (cf. Figure [TO]). 

A second reconstruction option is FLIC (Piecewise Linear Interface Calculation) - ||7l |33l 
and references therein. In this method, the interface is represented by a discontinuous chain of 
segments each obeying the definition [7]: 


xn = a. 


(52) 
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(b) PLIC 


Figure 10: The two main interface reconstruction approaches for the VOF method - SLIC and 
PLIC - are compared, (a) In the SLIC approach, the interface is reconstructed using 
line segments, which can he either parallel or perpendicular to the major flow axis, 
(h) In the PLIC approach, the interface is represented hy a discontinuous chain of 
segments. 


The parameter a can he determined using the condition of mass conservation. A central question 
for this scheme is the way, in which the normal is evaluated. The obvious choice, evaluating 
the normal the gradient of the volume fractions n = VF/|VF| approximated through finite 
differences as proposed in II223II does not perform well II1561I . which brought up a variety of 
alternative choices II179I . II 1561 presents a generalization of PLIC to adaptive moving grids, 11181 
to arbitrary unstructured grids. Note that, although based on piecewise linear functions, the PLIC 
approach cannot in general represent linear interfaces. II1881I proposes an appropriate extension. 
The definition of the orientation of the line segments can be enhanced by spline interpolation of 
the interface in every cell II142II . A special challenge is the treatment of more than two fluids, 
which connect at a triple point within a cell. In equilibrium, this point is connected to the contact 
angles of the particular fluids. In the instationary case, ll^ derives a relation to the normal 
vectors obtained within a PLIC setting. 

All piecewise linear reconstructions, i.e., both SLIC and PLIC, suffer from the unphysical 
creation of droplets disconnected from the actual surface by pure construction of the numerical 
method 

A similar approach, using a local height function, has been introduced in II127II and was for 
example utilized in II147II and references therein: The method operates on 3x3 blocks of cells. 
After the orientation of the surface has been determined based on the values of the volume fraction 
field of surrounding cells, the surface height is determined row- or column-wise. Then, instead of 
updating the volume fraction values, the height function is updated. A subsequent adjustment 
step of the corresponding volume fraction values is required. 

A very fruitful alternative is the coupling of VOF with the level-set approach as it will be 
described in Section]^- Coupled Level-Set Volume-of-Fluid (CLSVOF). Commencing with 
112051 . the idea has been propagated to combine the level-set approach - with its smooth interface 
- and the VOF approach - with its inherent mass conservation. Both the level-set function and 
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the volume fraction are advected simultaneously, yet individually. The volume fraction interface 
description is then the basis for a correction of the zero level-set fT). Such approaches, illustrated 
in Figurehave for example been explored in 112181 1^. A similar method is the reconstructed 
distance function (RDF) in 15411 . except that here the distance function is deduced from the 
piecewise linear reconstruction of the volume fractions. 



Figure 11: Coupled Level-Set Volume-of-Fluid: Both the level-set function and the volume 
fraction are advected simultaneously, yet individually. The volume fraction interface 
description is then the basis for a correction of the zero level-set (cf. Section |^. 

For some applications, simplified approaches that do not rely on a full interface construction 
can also lead to satisfactory results. For example, Q names WLIC (Weighted Linear Interface 
Calculation) 112221 1 150]| and THINC/SW (Tangent of Hyperbola Interface Capturing with Slope 
Weighting) II220I . WLIC is a simplification of FLIC. In a first step it treats each coordinate of the 
normal vector separately, thus computing fluxes along coordinate directions. In a second step 
these fluxes are averaged to account for slanted normal vectors. 

5.2 Curvature Computation 

The general definition of the curvature K in the VOF context is: 


However, being based on the heaviside function, the volume fraction function F is not smooth 
enough to be differentiated directly. Options can be found in either reconstruction of the interface 
(as described in the previous section), or some type of smoothing. 

An alternative formulation is: 


tc = Vn. 


(54) 
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Firstly, we will discuss methods for curvature computation that are based on interface recon¬ 
struction. In II127I . Kleefsman et al. suggest to compute the curvature from the height function h 
using the relation: 


K = 


dhjdy 


dy \ y/l -f [dhjdyY 


(55) 


The derivatives of the height function needed here are approximated using finite differences. 
Improved versions are among others offered hy lll43L who introduced a correction scheme based 
on local error estimation, and II1391I . whose scheme involves a collection of possible height 
functions and corresponding curvatures, which are then sampled using quality statistics. In 
the context of PLIC, II155II uses a local approximation of the curvature with polynomials. The 
polynomial coefficients are obtained from a least-squares fit. Raessi addresses the problem from 
a different perspective. In ||l84l . unit normals are advected, which are the principal constituent 
for both the reconstruction of the interface and the curvature calculation. 

Secondly, we will address the topic of smoothing. In Il54ll . three different approaches to 
curvature approximation are compared: (1) In convolved VOF (developed in ll21911 1. F is replaced 
by a smoothed version F, which has been mollified by a smoothing kernel K. K could for example 
be defined as K{r) = (1 — for r < 1. The danger here is that the smoothing can be so strong 
that the curvature effects are diminished. (2) An improvement with respect to error, but not 
with respect to convergence behavior, was found using RDF. (3) Superior to both method was 
however a height function approach by Sussman II204II . Another commonly used approach to 
circumventing interface reconstruction is to employ the continuum surface force (CSF) model 


as described in Section 3.2.2 Very recently, Baltussen has suggested the tensile force method 
fT5l . This is based on a substantially different surface tension model: the tensile force model 
introduced in ll214t . In this model, the surface tension force is determined as a sum of tensile 
forces, i.e., the force neighbouring cells exert on each interface element. This method does not 
require any curvature information. In |[60l a method the authors name CELESTE is presented. It 
revolves around a second-order Taylor expansion of F from any given cell to its neighbouring 
cell. It is constructed from an overdetermined system of equations, utilizing a least-squares fit to 
determine both the normals and the curvature. 

An interface capturing method which, contrary to VOE, provides a smooth interface is the 
level-set method, which will be subject of the next section. 


6 The Level-Set Method 

The level-set method is an interface capturing strategy. Level set methods 111731 |35l mitigate 
difficulties associated with a discontinuous volume-fraction function by using instead a smooth 
level set function (j ), separating areas filled with fluid A (0 < 0) from those filled with fluid B 
(0 > 0). The interface P"' is then given by: 


r'"' = {xea|(/)(x,t) = o}. 


(56) 
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In most cases, the level-set function <p is defined as a signed distance function II1721I : 


0(x,t) = ±min ||x —x*||, Vx G fl. (57) 

x*6r;"' 

At each point x and t, the level-set function (j) stores the shortest distance from the point to the 
current interface. Figure [T^ illustrates the concept by example of a circle. 




(a) Implicit circular inter- (b) Level-set field of a circular interface 

face 

Figure 12: A 2D domain containing a circular interface: In (a), an implicit definition of the 
interface is shown. Note that the elements of the mesh do not conform with the 
interface. The implicit description is realized through a signed distance level-set 
function outlined in (b). 

The method requires the definition of an initial level-set field from the original interface 
position: 


(/)(x,0) = (/)°(x) gD. (58) 

Furthermore, at inflow boundaries, the level-set field has to be defined throughout the entire 
time domain, thus assuring that the level-set function is well-defined at boundary points, where 
information is transported into the spatial domain. During the progress of the simulation, the 
initial (j) is then transported with the fluid velocity u(x,t) using the following transport equation: 

|^+u-V(/>=0 inD, fG[0,r]. (59) 

The value of the level-set function indicates, which density and viscosity values p, and p,- are 
associated with which element, or, along the interface, only part of an element. To illustrate the 
latter case. Figure [T^ shows a sample element, which is cut by the interface. The interface is 
approximated (e.g., linearly) through the intersection points with the element edges. Depending 
on the side of the interface where an integration point is located, the material property associated 
with it is chosen: 
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(a) Variables jump at the (approxi¬ 
mated) interface 



(b) Variables are smoothed across the 
interface 


Figure 13: A quadrilateral element with four nodes is intersected by the interface (red): The 
interface is approximated linearly (brown), (a) During integration, all Gauss points x 
are associated with either the properties of fluid 1 (pi,/ii) of of fluid 2 (P 2 ,P 2 )- (b) 
As an alternative, smooth transition of the properties can be generated throughout the 
element. Note that in both cases, an integration point might lie in the area between 
the interface as indicated by the level-set function and the linearized interface. One 
consequence of the approximation is then that this integration point will be associated 
with the wrong material domain. 




Pi,Pi, (/)(x,t)<0; 
P 2 ,P 2 , (/>(x,t)>0. 


(60) 


If standard discretization methods are used, this procedure may lead to stability problems. 
It is therefore suggested to create a smooth transition of the material properties through the 
element. This can for example be achieved by including a smoothed Heaviside function into 
the formulation 12061 . The price to pay is, however, a smeared out interface where it was once 
infinitesimally thin. 

Already in its original set-up, the level-set method attains the main advantage of the interface¬ 
capturing methods: the topological flexibility. Furthermore, as opposed to other interface 
capturing methods, the evaluation of geometrical quantities is straightforward. As (/> is a smooth 
function, it can be used directly to compute the interface curvature, normals, etc. What is left 
is to cope with the challenges that this approach contains: treatment of discontinuities across 
the interface, ensuring mass conservation, applying boundary conditions on the interface, and 
evaluating equations on the interface. 

Note that the volume-of-fluid method of the previous section can be interpreted using the same 
scheme, except that (j) = Heaviside. 


6.1 Mass conservation 

6.1.1 Reinitialization 


Since the velocity field entering the transport equation ( [59| l is usually non-uniform, it will over 
time degrade the signed-distance property of the level-set function. Even though the interface 
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is still represented correctly, this has severe influence on the quality of the computation of 
geometrical properties from the level-set function. Consequently, frequent reinitialization, i.e., 
restoring of the signed-distance character of 0, is recommended. 

Different approaches for reinitialization can he found in literature II115I . In direct reinitialization 
approaches, the closest distance to the interface is determined individually for each node. The 
approach is based on identifying the intersections between the interface and element edges. Out of 
the collection of these intersection points, the Euclidean distance to the nearest point is computed 
for each node. This value is stored as the new level-set value for the node, keeping the sign intact. 
Unfortunately, this approach can quickly lead to efficiency problems. The efficiency can be 
improved by restricting the reinitialization to a narrow band around the interface ll^ 1177112 1511 
(cf. Figure [14]). As an alternative, fast search tree methods can be employed to organize the 
distance evaluations ||149l . Fast marching methods compute the distances sequentially from low 
to high, each computation building on the previously computed distances. ||191I features timings 
for the different approaches. 




(a) Degraded signed distance function (b) Locally reinitialized signed dis¬ 

tance function 

Figure 14: To handle the degraded (as compared to Figure [T^ level-set function (a), a reinitial¬ 
ization restricted to a band around the interface is performed (b). 


6.1.2 Mass correction 


One very undesirable side-effect of the reinitialization procedure is the mass loss: The reinitial¬ 
ization does not preserve the location of the interface. This effect is enhanced through the choice 
of temporal and spatial discretization of the transport equation ( ^ i lll87t . Investigated solution 
approaches are manifold. Some approaches aim at the minimization of the displacement during 
reinitialization II 1031 1 163112061 . In |[52l an iterative mass correction is described, which provides 
global mass conservation. As mentioned in the Sections]^ andanother possibility lies in the 
combination with particle or volume-of-fluid methods. 


6.2 Obtaining sharp interfaces 

With the implicit definition through the level-set function, the discontinuity conditions along the 
interface wifi usually occur inside of elements, where, as described before, a smearing of the bulk 
properties is recommended. One possibility that allows to account for a sharp interface within the 
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element is the extended finite element method (XFEM). XFEM is then often connected to the 
other improvement option: adaptive mesh refinement. 

6.2.1 XFEM - the extended finite element method 

In the finite element method, the unknown function is usually interpolated with polynomial basis 
functions that are continuous within the elements and continuous across element faces. This 
approach is very much suitable for smooth unknown functions, but leads to problems if jumps 
occur. The XFEM circumvents this deficiency by local and well-targeted enrichment of the basis 
function space. The origin of XFEM lies in fracture mechanics I161II . Overviews of the 
method can be found in In XFEM, a sample function g is approximated on a fixed mesh 

as: 


/(x,0 + (61) 

iei iei* 

' -V-" '-V-' 

standard FE enrichment 

Ni{x) and N*{x) are the standard finite element shape functions for node i with the nodal un¬ 
known gi. The set I is the set of all nodes in the domain, whereas 1* contains all enrichment nodes. 
They are enriched using the enrichment functions v^(x,t), which are connected to additional 
XFEM unknowns a;. Refer to Figure [T^ for a sample mesh. 


. I 

□ 7 * 

blending elements 
^ enriched elements 
— interface 

Figure 15: Sample domain with an interface. The three different element types of the XFEM: 
regular, enriched, and blending elements, are indicated. 

An enrichment function that is typically chosen for strong discontinuities is the sign-enrichment: 



Wsign{^,t) = sign{^{x,t)) 


'-1, 

(/>(x,t) 

<0, 


0, 

(j){x,l) 

= 0, 

(62) 


(j){x,t) 

>0. 



In the fluid flow problems considered in this paper, the scenario is as follows: Across the 
interface, a jump in density and viscosity, a kink in the velocity, as well as both a jump and kink 
in the pressure need to be considered (cf. Figure]^. 

Chessa and Belytschko IIM1I3911 use a kink enrichment for only the velocity approximation. 
Depending on whether surface tension is considered, a kink or jump enrichment of the pressure 
is added 117^1159111851122511 . Koike 1113111 enriches both spaces with a jump. Even so it seems 
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like a natural enrichment to account for the kink in the velocity field hy appropriately enriching 
the basis, B9l states that this does not lead to an improved overall quality of the solution. 111921 
even reports stability problems if this is included. As a consequence, several authors apply an 
enrichment for only the pressure IITTI 1^1951 . 

It is to be noted that the enrichment also entails special quadrature procedures. In order to be 
able to invoke standard Gauss rules, the interface cells are subdivided into triangles as indicated 
in Figure El|79lmu. 



(a) Original interface 



face 



(c) Sub-cells and 
Gauss points 


Figure 16: Linearization of a curved interface in a quadrilateral reference element by decomposing 
the quadrilateral into linear triangular sub-cells, where finally fhe infegrafion poinfs 
are placed. 


6.2.2 Adaptive mesh refinement 


An alfernafive or an addition fo XFEM can be found in adapfive mesh refinemenl (AMR). AMR 
can be used to increase the solution quality in the vicinity of the interface: the most interesting 
part of the computational domain. While global mesh refinement usually leads to prohibitively 
large mesh sizes, AMR provides an alternative at reasonable computational cost l[T4lll801l . In 
lIT/llSOl . an adaptive refinement approach tailored to XFEM is introduced. Flere, the level-set 
functions serves as an indicator for the refinement area. The interface elements are successively 
subdivided up to a desired level. Eigure[^ depicts the general refinement procedure for a simple 
setting. Flere, two levels of refinement are applied. In a first step, the elements of the initial 


mesh. Figure 17(a) which are in contact with the interface, are subdivided once (cf. Figure 
|17(b) I. Those elements from the first refinement, which are cut by the interface are refined a 
second time, see Figure 17(c)| Subsequently, further elements are subdivided step by step such 
that neighbouring elements differ by at most one level of refinement (cf. Figure 17(d)| ). This 
assures a smooth variation of the element sizes. The level-set indicator field is defined on a fine 
background mesh with a mesh density corresponding to the highest refinement level (cf. Figure 
|17(e) I. Thereby, no interpolation of the level-set values at nodes created during refinement is 
required and the accuracy of the interface does not degrade. 


6.3 Computing geometrical quantities from the interface 

Since the level-set function (p is smooth, it can be used to compute geometric entities. As an 
alternative, the (linear) approximation of the level-set function within the elements could be 


29 
























(a) Initial mesh 



(b) One level of refinement in 
interface elements 



(c) Two levels of refinement 
in interface elements 
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(d) Additional refinement in (e) Level-set indicator field 

neighbouring elements (interface) defined on a back¬ 

ground mesh 


Figure 17: Adaptive mesh refinement with respect to the interface (bold red line). 


employed. In terms of convergence order, ll^ 1^ suggest the second option. 

6.3.1 Normal vectors 

The normal vector is defined as: 


n = 


V0(x,f) 




(63) 


Due fo fhe signed disfance properfy of 0, fhe denominator in ( [6^ could essenfially be dropped. 
Flowever, since during the time-stepping, the level-set function is only an approximation of a 
signed distance function (cf. Section [67T]), it is advisable to keep the normalization. 


6.3.2 Curvature 

The curvature can be computed from the level-set function as: 


V(/>(x,t) 


V(/>(x,t) 


pi/ir 


(64) 


This formulation is restricted to approaches, where the level-set function is not approximated 
linearly. An example is given in Il70l . where a quadratic approximation in combination with a 
filtering technique is employed. However, even when theoretically possible, Marchandise points 
out in II149I that it is in general not advisable to use the thereby computed curvature to evaluate 
the surface tension. One popular option is to resort to the Laplace-Beltrami technique already 
described in Section [3. 2. II 

Particularly in combination with the XFEM, level-set provides very sharp interface descrip¬ 
tions. A contrasting approach is the phase-field approach described in the next section, which 
deliberately introduces diffuse interfaces. 
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Figure 18: Illustration of the phase-field variable: ^{x,t) takes on one value (e.g., 1 in this case) 
in one phase and another value (in this case —1) in the second phase. The transition is 
rapid, but not sharp. We speak of a diffuse interface with interface width 5/. 


7 Phase-Field Method 

One of the oldest numerical methods in multi-phase problems is the phase-field method. Ac¬ 
cording to f/Ol . as early as 1873, the work of Gibbs on thermodynamics already served as a 
foundation |f89l . 

The main difference compared to the previously introduced interface-capturing methods is 
that the phase-field mefhod works with diffuse interfaces — i.e., the transition layer between the 
phases has a finite size. There is no tracking mechanism for the interface, but the phase state is 
included implicitly in the governing equations. The interface is associated with a smooth, but 
highly localized variation of the so-called phase-field variable (p. In two-phase problems, 0 is a 
scalar value; if more phases are present, it can become vector-valued. <p is assigned a distinct 
value, e.g., —1, in phase A and another distinct value, e.g., 1, in the other phase B. The interface 
could then be assumed to be located at (/> = 0 ll30l ; however, the phase-held equations never 
require the knowledge of the exact interface location. ^ will vary with space and time. Typically, 
we observe small variations in the bulk phases and rapid variations close to the interface fTOl . 

Comparing Figure [TS] and Figure [T^ gives the notion, that phase-held resembles level-set in 
many ways. Indeed, they have been applied for identical applications, however, there are some 
fundamental differences that should be considered in the choice of application. Major differences 
are: (a) The phase-held approach avoids setting continuum boundary conditions at the interface, 
which would usually be required to match the individual solutions of the bulk phases. Note also 
that one of these boundary conditions usually dehnes the advection velocity of the interface. 
What is particularly difficult about these boundary conditions is that they have to be set at a 
boundary whose shape and position is also part of the solution. Instead, these boundary conditions 
are substituted by the inclusion of the phase-held variable into the system. Consequently, the 
phase-held method is particularly useful in cases, where the continuum equations for a sharp 
interface do not (yet) exist. Sharp-interface descriptions can be interpreted as a phase-held 
approach with inhnitely small interface width. However, it is not considered to be very useful to 
use phase-held in such a sense llTOl . (b) In the phase-held method, the interface is not advected 
specihcally as would be the case for level-set; moreover, the interface location is never computed 
explicitly llTOll^ . 
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Phase-field can be seen as a computational method, in which the presence of a non-sharp 
interface dampens instabilities, or it can be seen as a physically motivated approach, where the 
governing equations do not impose sharp interface conditions ifTOll . 

The derivation of the governing equations for a system described by a phase-field approach is 
strongly tied to the notion of thermodynamic equilibrium. In engineering applications, we usually 
deal with open or closed systems; and hardly with isolated systems. In open or closed systems, 
the thermodynamic equilibrium is not defined through the state of maximal entropy, as it would 
be for isolated systems, but through the extrema of some other thermodynamic state variable. In 
principle, we distinguish between four variations of the equilibrium by thermodynamic forces: 
mechanical, thermal, material, and chemical. 111461 

• Mechanical'. Two systems under different pressures are inclined to reach mechanical 
equilibrium: mechanical work is performed. In such a case, the pressure is interpreted as a 
thermodynamic potential. 

• Thermal: Two systems with different temperature are inclined to reach thermal equilibrium: 
heat flow is induced. In such a case, the temperature is interpreted as a thermodynamic 
potential. 

• Material: A system containing either different material states of one material, e.g., solid 
and liquid phase, or two chemically inert materials, will tend to the material equilibrium. 
In such a case, the chemical potential is interpreted as the thermodynamic potential. 

• Chemical: A system containing substances, which are chemically reactive, will develop 
a specific equilibrium between educts and products. Again, the chemical potential is 
interpreted as the thermodynamic potential. 

The state of equilibrium is characterized through the minimum of the thermodynamic potential 
P. The condition for equilibrium is defined as dP = 0. 

So far, we have stated that in order to describe a thermodynamic equilibrium of any kind, 
the definition of a suitable thermodynamic potential is essential. In realistic problems, several 
potentials interact. Therefore, the need for finding the extremum of one specific potential is 
now replaced by minimizing the free energy functional of the system. The free energy is the 
portion of the energy that is available to perform thermodynamic work. After performing such 
work, it is irreversibly lost fTOl . Depending on the applications, different definitions of a free- 
energy functional can be used: the Gibbs free energy formulation is for example applicable to 
systems where the temperature and the pressure remain constant, whereas the Helmholtz free 
energy formulation can be employed for situations where the temperature and the volume remain 
constant. In general we can write f/Ol : 



(65) 


where P refers to the appropriate thermodynamic potential with P signifying the respective density 
function. Xi,A 2 ,...,A„ refer to the governing variables of the system, i.e., the independent 
state variable. Depending on the system, these could for example be temperature, pressure, the 
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Figure 19: Double-well potential: A potential depending on the phase-field variable 0 is defined, 
which fulfills fhe basic criteria of symmetry and two minima. 


concentration of the phases, etc.. The free-energy functional together with the equations of motion 
leads to a full description of the system in the sense that the value of all state variables can be 
determined from the free-energy via extremal principles if we assume thermodynamic equilibrium. 
In detail, this means that through differentiation with respect to the individual variables, equations 
for all those variables can be obtained (cf. ifTOl and references Ill71[|57ll99lll81ll therein). In the 
phase-field method, an additional, in a sense artificial, dependence on the phase-field variable is 
generated. Equation ( [bb) ) now reads iTTOll : 

P(Xi,X2,...,X„,(/.)= f P{XuX2,...,X„,^)dV, (bb) 

Jv 

where (p is the phase-field variable. 

In order to include (j) into the formulation, a contribution of (j) to the density function P has 
to be dehned: the potential F(0). V((j)) shall be constructed in such a way that it features two 
minima - one for each of the two phases. This ensures that exactly these two phases (e.g., liquid 
and solid) are the most stable states. Furthermore, the potential needs to be symmetric. If we 
expand V((j)) in a Taylor series, and break off after the first two terms, we will find the definition 
|7Ql: 


= (67) 

Flere, it is assumed that the two minima have an energy level of 0. In order to be able to 
also represent inbomogeneous states, tbe potential also bas to depend on tbe gradient of (j). Tbe 
lowest order terms invariant under both rotation and translation are either or The 

contribution of (j) to the energy can then be expressed as: 

p = --- + ^iy(t>f+v{(t>), (b8) 

with £ as a constant. This formulation is termed the Cahn-Hilliard potential. 

Again using variational principles, evolution equations can be derived for the phase field. 

In non-equilibrium slate, two important equations that can be developed from the Gibbs 
free-energy functional are the Cahn-Hilliard and the Allen-Cahn equation. The latter is a reaction- 
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diffusion equation and in general not of importance for flow problems. The Cahn-Hilliard 
equation reads Il30l 


| = V.(Mcv(|f-<^Ac)), (69) 

In the above equation, the unknown quantity is the concentration field C of one of the two 
quantities present in the simulation (e.g., fluid A and fluid B or fluid and solid). C can be 
interpreted as equivalent to 0. /(C,...) denotes a free energy. It depends on the concentration, 
but also on other independent state variables such as possibly the temperature. Me is the mobility 
related to the solute diffusion coefficient and Ec is a constant related to the length scale of 
the diffuse interface. Note that the Cahn-Hilliard equation contains a fourth derivative of the 
concentration field; a challenge to many discretization methods. 

The most important areas of application of the phase-field method are solidification processes, 
especially including phase segregation in binary alloys, and dentritic growth llC71l9^[T98lll991 
120311 . It has also been applied to other scenarios, such as spinodal decomposition - the unmixing 
of a mixture of liquids or solids - ifT^ . surfactant transport EUTTl, or even topology optimization 


In 11I40I . isogeometric analysis (cf. Section 8.4.2) is utilized to analyze liquid-vapor phase 
transition phenomena modelled through the Navier-Stokes-Korteweg equations, which can be 
categorized as a phase-field model; in ll9ll the same approach is used for the solution of the 
Cahn-Hilliard equations in the context of spinodal decomposition. For both applications, the iso¬ 
geometric analysis approach profits from the arhitrarily high continuity of the finite-element basis 
functions across element interfaces. For the first time, this allows for a standard use of the finite 
element method not augmented by mixed methods or discontinuous or continuous/discontinuous 
Galerkin methods ||9T1. 

The discussion of the phase-field method concludes the topic of interface-capturing approaches. 
In the next section, interface tracking will be introduced. 


8 The Interface-Tracking Approach 

In interface tracking, the boundary or interface is explicitly resolved by the computational mesh. 
The necessity of accurately following the deforming boundary of the domain requires some 
degree of Lagrangian description to be present in the formulation, at least in the vicinity of the 
boundary. Since the use of a Lagrangian viewpoint is unavoidable, some researchers considered 
fully Lagrangian formulation for the entire domain 111041116811183111061116811 . The simplicity of 
such solutions is offset by several difficulties. The flow field in the interior of the domain may 
exhibit a large amount of circulation or shear, which are often not directly related to the motion 
of the boundary. The fully Lagrangian approach then results in unnecessary mesh distortion and 
invalid meshes in such interior regions, even for moderate or null displacements of the boundary. 
However, internal circulation and shear are handled without difficulty by an Eulerian description. 
This provided a strong motivation for the development of methods capable of combining the 
Lagrangian and Eulerian approaches in the same domain. The ALE (Arbitrary Lagrangian- 
Eulerian) formulation was initially stated in the finite difference context by Hirt 11107L and later 
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adopted also in the finite element community—see llhTll25l 11121 IllOH and references contained 
therein. The basic idea of ALE is to work with a deforming mesh, as in the Lagrangian method, 
hut to decouple the mesh deformation from the displacement of the fluid particles - or in other 
words, the mesh velocity is no longer equal to the fluid velocity, hut can he chosen arbitrarily. 
The usual choice in ALE is to track boundary and interface nodes in a Lagrangian manner 
(or a manner that comes very close), while all other mesh nodes are adjusted solely with the 
intention of preserving mesh quality as best as possible. The ALE approach requires that the 
mesh velocity enters the momentum equation ([T]l explicitly, thus modifying the original flow 
equations, specifically the convection term. Eurthermore, this equation now needs to be written 
over a reference domain that may or may not coincide with either material or spatial domains. 
This added complexity is one of the main drawbacks of AEE. Nevertheless, it is used actively in 
the area of fluid flow computation, e.g., in ll^l82l [851121711 . The AEE philosophy found another 
expression in the family of space-time finite element methods, which are described further below. 

8.1 Deforming spatial domain space-time finite element methods 

The classical choice of space and time discretization within the finite element community is to 
use finite elements (EE) in the spatial domain, but finite differences (ED) in the time direction. 
One alternative are the so called space-time finite element methods, which use finite element 
approximations for both space and time. Eeaving the spatial discretization unaltered as compared 
to the classical approach, an additional coordinate direction - the time t - is introduced for both 
the computational domain as a whole, but also for each finite element. As a consequence, all 
volume integrals in Equation ( [T^ are now integrals not only over the spatial domain D, but over 
the space-time domain Q of dimension nsd + 1. The shape functions A;, which interpolate the 
unknown functions, are equipped with an additional dependency on the time: A(x,t). 

Although it is in principle possible to employ the space-time discretization with an interpolation, 
which is continuous-in-time, the modern space-time approach usually relies on discontinuous in 
time discretization in the spirit of the Discontinuous Galerkin method lIMllllllllOlll . In order to 
avoid unmanageable number of degrees of freedom that must be solved for at any given time, 
the space-time approach is typically applied to subsets of the temporal domain called space-time 
slabs (cf. Eigurej^, which are roughly comparable to time steps in the EE/FD approach. Due to 
the discontinuous in time interpolation, each node stores two values of every unknown per point 
in time. In order to enforce a relation between these two values, the continuity of the solution in 
temporal direction is imposed weakly by adding a jump term to the variational form ( [T^ . 

The space-time method offers the unique opportunity to write the variational form directly over 
a deforming domain. This concept was pioneered by Jamet and Bonnerot II121[|^ 177111481112011 . 
It was later extended concurrently by Tezduyar et al. II207II in the form of the Deformable- 
Spatial-Domain/Stabilized Space-Time (DSD/SST) method, also applied in Il23lll601l . and by 
Hansbo II 1001 with a particular mesh moving scheme conforming to the characteristic streamline 
diffusion ideology. Note that in the deforming domain space-time approach, there is no need to 
explicitly include the mesh velocity into the variational form, as it would be mandatory for AEE. 

In most space-time implementations to date, the meshes for the space-time slabs are simply 
extruded in the temporal direction from a spatial mesh, resulting in reference element domains 
that are always Cartesian products of spatial and temporal reference domains. Such an approach is 
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Figure 20: Illustration of two successive space-time slabs From the lower time level to the 
upper time level, the mesh can deform. In blue, one space-time element is exemplified. 
It is triangular in space and forms a prism in space-time. 


best described as semi-structured (unstructured in space, structured in time) and does not leave the 
option of increasing temporal refinement in portions of the domain. In such a case, the space-time 
slab exactly corresponds to a time step of a semi-discrete procedure. In fact, many stencils of 
the FE/FD methods may be re-derived by using the semi-structured space-time approach with 
appropriate weighting and interpolation functions. 

Note that the applicability of the space-time approach is of course not limited to the de¬ 
forming domains in an interface-tracking context. Extensions of the XEEM concept to space- 
time Il40lll751ll51l . the monolithic space-time fluid-structure-interaction solvers II109II . or shape 
optimization ||6^I176II are just a few examples. 

8.2 The mesh deformation 

In interface tracking, the mesh is continuously adapted to the flow. Procedures for mesh deforma¬ 
tion are described in the following, categorized into boundary deformation and displacement of 
the inner nodes. 


8.2.1 Displacement of the boundaries 


Section 3.1.2| describes the standard choice for the motion of the domain interface/boundary in 
the interface tracking context: a full Lagrangian deformation as specified by Equafion ( [2T] ). This 
franslafes fo a deformation, where fhe mesh velocify v is chosen exactly equal to the fluid velocity 


v(x)=u(x). (70) 

However, this is not the only possible choice. Physically, the deformation of the inter¬ 
face/boundary mimics a no-penetration boundary condition, meaning that no fluid is allowed to 
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Figure 21: Possible displacement directions satisfying the kinematic boundary condition. 


cross the interface - this condition will always be fulfilled if the interface moves in accordance 
with the fluid as in Equation ( f/Ol ). However, it will also be fulfilled by any choice where the mass 
flux through the boundary, m, is 0. The mass flux m through the interface can be computed as 

m- 


m= / p(u —v)ndx = 0. (71) 

7r"" 

Consequently, all choices for v that comply with the kinematic boundary condition 


v(x) • n(x) = u(x) • n(x) 


(72) 


are valid (cf. Figure [2T]). It can be deduced from Equation ( f/I] ), that tangential node movement 
remains irrelevant to the interface/boundary shape, but may significantly influence the mesh 
quality. Particularly applications with large tangential velocities - such as die swell in extrusion 
processes or rising droplets - profit tremendously if the tangential velocity component is modified 
or even suppressed when computing the interface/boundary velocity v. 

Note, however, that the assumptions made by the kinematic boundary condition are only valid 
strictly for the analytical geometry of the interface/boundary. This is particularly due to the 
computation of the normal vector, which is usually ambiguous at the finite element nodes (usual 
finite element discretizations employ shape functions, which offer only C** continuity across 
element boundaries). Still, the finite element nodes are exactly the points of the interface/boundary 
where the mesh velocity needs to be defined. It is therefore inherent to the nature of the 
discretization that Equation (72 1 can only be fulfilled approximately by any choice other than 
Equation ( |70| ). II202II 

Nevertheless, Behr ll^ indicates several possibilities for boundary displacement resulting 
from Equation (|7^: displacement with the normal velocity component with v = (u • n)n and 


displacement only in a specific coordinate direction ( e.g., y-direction), i.e., v = Note here 

that a direct projection onto the unit vector Cy in the coordinate direction perpendicular to the 
main flow direction is not possible, as this would result in a displacement in the correct direction, 
but not with the correct magnitude to satisfy the kinematic boundary condition. 
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In cases with very extensive deformations, even this method can lead to strongly distorted 
elements. As a remedy, it is possible to allow the nodes to slip along the tangential direction. 
The amount of tangential movement is then not connected to the fluid velocity, hut determined 
according to the same method responsible for the inner nodes, which will be described in the next 
section. 

8.2.2 Retaining mesh quality: mesh update for inner nodes 

As described in the previous section, both the ALE approach and the deforming domain space- 
time approach prescribe only the mesh displacement on (or possibly in the vicinity of) the 
boundary/interface. All other nodes have to be adapted to this displacement in a way that retains 
good mesh quality. Wall 1121611 identifies two possible categories for the treatment of the inner 
nodes: methods that rely on explicit functions, which have been defined a priori, and methods 
that require the solution of their own system of equations. 

Examples of the first category can be found in II182[|126]| . where the boundary deformation is 
distributed smoothly throughout the entire domain using an interpolation kernel. 

In the second category, probably the most popular option is the Elastic Mesh Update Method 
(EMUM) HI221 and comparable approaches (e.g., Ill521ll66l f. In this method, the computational 
mesh is treated as an elastic body reacting to the boundary deformation applied to it. Eor this 
purpose, the linear elasticity equation is solved for the mesh displacement 1), which relates to the 
mesh velocity v as 1) = \At: 


^ ■ ®^mesh 0 , 

O^mesh(l^) —(h'^mesh(l^)) I T 2^ 
£mesh(u) = ^ (Vl) + (Vl))^) . 


(73) 

(74) 

(75) 


In structural mechanics, X,„esh and Hmesh are the Lame-parameters that are needed to define the 
stiffness tensor for isotropic and linear elastic materials. In the context of the elastic mesh update, 
these parameters are prescribed for each element in order to control its stiffness. This is used to 
increase the stiffness of the smaller elements compared to the larger elements in order to allow 
for larger deformation before the mesh becomes invalid. The choice of Lame-parameters can 
differ significantly between the approaches. 

As boundary conditions, the displacement of the boundary/interface nodes, which is prescribed 
due to the rules given in Section 8.2. 1[ is imposed: 


t) = 1) on rt)(t). 


(76) 


Another possible boundary condition, which, in the context of mesh deformation, is usually 
used in addition to the condition above, would be: 


II ■ ^mesh — hmesh OH ^hmesh (0 • 


(77) 
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(a) Undeformed mesh (b) Elastically deformed mesh 

Figure 22: Elastic Mesh Update Method: An initially circular interface (a) is subject to a non- 
uniform motion, (h) shows the ideal elastic response of the entire mesh with respect to 
the interface movement. 


This condition lets the nodes slip along certain boundaries. ru(t) and no longer 

distinct subsets of r'”'(t), but can be combined for different degrees of freedom. The calculated 
displacement l) can now be used to update the nodal coordinates in all phases: 


x«+i^x" + i). (78) 

As a possible alternative, a simple Laplace equation can be solved for the displacement instead 
of Equation ( |7^ , which diffuses the displacement throughout the domain, but gives less control 
over the individual elements. 

Also in the second category, but from a completely different view point, a method where a 
functional associated with the mesh distortion is minimized is proposed in II141II . The approach 
presented in HI3011 is based on a series of optimization steps with regard to the mesh quality. 
Loubere et al. 1 14411 apply a method, where all nodes are displaced in a fully Lagrangian fashion, 
but each time-step contains a reconnection step, where a new connectivity is established based on 
the existing nodal positions. 

As a final option, remeshing at every time-step can be applied, as for example in 111671 in 
the context of the particle finite element method (PEEM). Strictly speaking PEEM is a particle 
method, with marker particles carrying the information of the physical properties. What relates 
this method to the interface-tracking approach however is that the particles are connected by 
a Lagrangian mesh, where the marker particles simultaneously act as finite element nodes, on 
which the flow equations are evaluated. 

8.3 Normal vector and curvature approximation 

In principle, the explicit description of the interface/boundary allows the immediate evaluation of 
normals and curvature within the elements, as long as the shape functions can be differentiated 
sufficiently often. If the values for the normal vector or the curvature are needed directly at 
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element nodes, they can be obtained by a weighted average of the values within neighbouring 
elements. One example for defining this average can be found in ITTHI . 

In cases, where the geometric interpolation cannot provide the necessary geometric information 
- i.e., mainly in the case of linear interpolation -, one may either resort to the Laplace-Beltrami 
technique (cf Section 3.2.11 or introduce an additional geometrical entity from which the desired 
values can be obtained. Here, Mier-Torrecilla et al. lll58t propose the use of the osculating circle 
to the curve at each grid point. In the implementation, the osculating circle is identified as the 
circle with shoulder points at the respective node and its two neighbours. Using the radius r of 
the osculating circle, the term Jcn can be evaluated as y^. The use of spline interpolations for 
this purpose is proposed in 1821 with the use of piecewise defined cubic splines, and in |[6^ with 
the use of a single NURBS curve or surface to represent the boundary. 

The utilization of splines to obtain geometric quantities in interface tracking naturally leads 
to methods, which integrate Computer-Aided-Design (CAD) information into the finite element 
analysis. 


8.4 Recent advancement in CAD-integration for the finite element method 

Geometries for engineering applications are generated using Computer-Aided-Design (CAD) 
systems: the CAD model is what we assume to be the real geometry. As of now, all major CAD 
systems share one common basis for geometry representation: Non-Uniform Rational B-Splines 
(NURBS). They provide a a common standard to exchange geometry information. Similar to the 
close relation between level-set methods and XFEM, an interconnection between CAD-related 
numerical methods and interface tracking can be identified. (Nofe however thaf in a few cases, 
the level-set method has also been combined with spline descriptions, as for example in ll27ll8Ti .') 

In the classical discretization methods, the exact NURBS geometry is lost as a finite element 
mesh, which is in general only an approximation of the geometry, is generated. In this section, 
we will explain the structure of NURBS and then introduce two methods that make use of the 
NURBS format in order to integrate the exact geometry into the finite element method, thus 
avoiding the approximation character of the finite element mesh. 


8.4.1 Geometry representation using splines 

We will start out with the most basic NURBS structure: the NURBS curve. From the curve defi¬ 
nition, it will later be straightforward to derive the definitions of higher-dimensional objects. This 
section is based on ||189[I178]| . NURBS are an example of a parametric geometry representation. 
In this type of representation, a local parameter (e.g., 6) is defined along the curve. Usually, we 
use 6 = 0 at one end of the curve and 6 = 1 at the other end. The global coordinates of the curve 
{x,y,z) are then each defined as functions of 6: 


x = x(d), y=y{e), z = z{Q). (79) 

Important specifications during the historical development of NURBS were to generate a curve 
(1) whose smoothness is completely under user control. 
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( 2 ) which occupies a predefined space, 

(3) which allows for local shape control, and 

(4) which is able to represent both free-forms and analytical shapes. 

In the NURBS context, the space occupied by the curve is defined using control points, which 
are connected to a control polygon. It is a property of the NURBS that it will always remain 
within the convex hull of the control points (Requirement (2)). By modifying the position of the 
individual control points, the shape of the curve will also be modified, as the control points guide 
the curve. The NURBS definition allows giving certain control points a larger influence on the 
curve compared to others, as each control point is assigned a weight. Another important feature 
of a NURBS is that each control point is only responsible for the shape of a minor portion of the 
curve (Requirement (3)). In the NURBS definition, this effect is realized through the B-spline 
basis functions M, the building blocks of the NURBS basis which will be introduced below. 
The B-spline basis functions are non-zero only in a sub-domain of the parameter space. Each 
control point i is associated with exactly one basis function M,-p of polynomial degree p. It is the 
non-zero domain of this basis function that indicates which portion of the curve is influenced by 
the specific control point. This indication is given in terms of the local parameter 6. A typical 
basis function for a NURBS is depicted in Figure]^ The control point associated with this basis 
function will only influence curve points associated with 0i < 6 < 64 . The values of 6, where 
the influence of a basis function begins and ends, are referred to as knots. 



Figure 23: Example of a B-spline basis function with polynomial degree 2. Mi 2 is non-zero 
between 61 and 64 . 

The curve definition C of a NURBS curve is then 


c(e) = £Ri,p(e)Pi, ( 80 ) 

i=l 

where denote the NURBS basis functions and P, are the coordinates of the control point i - 
which are always given in the global coordinate system (x,y,z). n indicates the total number of 
control points. As we can see, C is a mapping from the parametric coordinate 6 to the global 
coordinate system (x,y,z). Notice that a NURBS curve has only one local parameter, no matter if 
it exists in or whether it is embedded into the or the M^. In this respect, NURBS are an 
example of an immersion. 

The NURBS basis functions R emanate from the B-spline basis functions M. These in turn are 
usually defined recursively (cf. Figure [24]), starting from the constant basis functions, using the 
Cox-deBoor relation: 
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Figure 24: Recursive definition of the B-spline basis functions: For each degree p, the total 
number of basis functions reduces by one. Each basis function of the higher degrees is 
built-up of p -f 1 constant basis functions, meaning that it is non-zero exactly between 
6i and 0,+p+i. In any interval 6i < 6 < 6i+\, exactly p -|- 1 basis functions will be 
non-zero. 




1 , ei<o<Oi+y, 

0 , otherwise , 


(81) 




d-di 


(82) 


R is then defined as: 


Ri.pW 


Mi^p{e)wi 


(83) 


The rational structure of the basis is responsible for the fact that NURBS are able to represent 
all analytical shapes, including conic sections (Requirement (4)). The values 6,, Qi^p, etc. are the 
knot values, i.e., the indication of the influence domain of the basis functions. The knot values 
are collected in the so-called knot vector: 


0 = [^ 1 , ^ 2 , ^ 3 ) ^ 4 , • ■ •] 


Mi,o^O Ms.o^O 


(84) 
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For the proper definition of the basis, it is essential that the knots are in a non-decreasing 
sequence. However, it is specifically possible to repeat knot values. With each repetition, the 
continuity of the basis is decreased. With regard to the curve, there are then two measures 
that influence the continuity: the continuity of the basis and the topology of the control points. 
If control points are placed on top of each other, the curve continuity is decreased, e.g., to 
represent sharp comers. If control points are aligned, the continuity can be increased. Both 
effects are a direct consequence of property that the NURBS curve is always contained within 
the convex hull of its control polygon. These two measures give complete control over the curve 
continuity (Requirement (1)). In general, the basis functions, and therefore the curve, are infinitely 
differentiable at all points, and p — 1 times differentiable at the knots. 

The derivatives and the curvature of NURBS curves can be calculated analytically. From 
differential geometry, the curvature of a parametric curve C(0) can be calculated with ||94l: 


K(e) 


c'(e)xc"(e)i 

“lew 


(85) 


C'(6) and C"(6) denote the first and second derivative of the curve with respect to parameter 6. 
For the sake of brevity, Mi p{6) is written as and as in the following equations. 


c'{e) 

c"{e) 


LiM'i pWiPiY^iMi^pWi - LiMppWiPiZiM- pWi 
{LiMi,pWif ’ 

(LiMi^pWif 


21; M'i pWi M'f pWiPi I; Mi^pW, 

iLiMi,pWif 


2 Li M'i,pWi Mi^pWiPi Li M'i pW, 

{LiMupWif 


( 86 ) 


(87) 


M'- p and M-^ are the first and second derivative of the basis functions with respect to the parameter 

e.’ 

Figurej^gives an example of a NURBS curve representing a circle. A closed curve is obtained 
by overlapping the first and the last control point. Higher continuity at the connection point of the 
curve can be reached if p control points at the ends are overlapped (P„_p+i = Pi,..., P„ = Pp). 

If surfaces (or volumes) are to be described through NURBS, a second (or third) local parameter 
is introduced. The basis functions for the higher-dimensional objects are computed as a Cartesian 
product of the univariate basis functions M/,p(0) and Ty^(i) defined in Equation ( [^ . For each 
local coordinate direction, an individual knot vector is required. Note that it is possible to choose 
different polynomial degrees of the basis for each local coordinate. The NURBS surface S is 
described as 


S(e,l) = ££<l 0 ,l)P/,;, 


!= 1;=1 


( 88 ) 
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Pa=(-l,1.0) P5=(0,l,0) P4=(l,l,0) 



P8=(-l,-l,0) P,=P9=(0,-1,0) P2=(lrl,0) 

Figure 25: Example of a NURBS curve: with 9 control points and a quadratic basis, it 
is possible to represent a full circle. The corresponding knot vector is 0 = 

[0 001/4 1/41/21/23/43/41 1 1]. The weights have to be set to wi = W 3 = W 5 = 

/2 

W 7 = W 9 = 1 and W 2 = W 4 = wg = W 8 = 


with 


R//(e,i) = 






(89) 


8.4.2 Isogeometric analysis 


Recall that the finite element method invokes the isoparametric concept, as described in II113II 
or many other works on finite elements. The isoparametric concept signifies that the element 
nodes are interpolated with the same shape functions as the unknown function, i.e., the solution 
of the partial differential equation. This approach instantly presents a mapping from reference 
coordinates to global coordinates, thus offering a simple means for evaluating all components of 
the transformation theorem necessary to compute the integrals of the weak form on a reference 
element instead of the global elements. In comparison to an integral evaluation in global 
coordinates, the approach using a reference element makes easy use of the local support of 
the shape functions and furthermore avoids the more tedious and computationally inefficient 
definition of the shape functions in global coordinates. 

Cottrell, Bazilevs and Hughes used the exact same underlying approach when devising the 
concept of isogeometric analysis Illl4[l50ll20ll . In the spirit of the isoparametric concept, it 
utilizes the NURBS basis functions in order to represent both the geometry and the unknown 


solution. For the geometry, the already presented definition ( [ 8 ()| ) or ( [ 88 ] ) is employed. The 
unknown function u is then approximated in exactly the same fashion as: 
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Pi P 2 P 3 P 4 

•-•-•-• -► ;c 

I-► e 

Figure 26: Function representation in IGA: Consider a quadratic NURBS curve whose control 
points are aligned in such a way that it represents a line (lower part of the picture). 
The local parameter along the curve is 6 , and it is embedded into with coordinate x. 
On this curve, a function u is approximated (u^). The function is represented through 
the control variables ui through U 4 which are interpolated using the NURBS basis 
function Ri p. 


u^ = Z^i,p{e)ui. (90) 

1=1 

The values m, are the control variables. These are the unknown values solved for in the finite 
element code. Note that just like the control points for the geometry, the values of the control 
variables are in general not coinciding with the solution; see Figure [26| 

Not only does this idea essentially provide the capability to perform finite element analysis 
directly on CAD models, but it also profits from superior approximation properties of NURBS as 
compared to the polynomial shape functions used in the standard finite element method. 

8.4.3 The NURBS-enhanced finite element method 

One major challenge still remaining for the isogeometric analysis described in the previous section 
is the generation of closed volume splines describing complex geometries. In the traditional 
CAD/CAM systems, complex geometries are usually described as a large number of individual 
surface patches, connected neither geometrically nor parametrically. The reason for this is that 
traditionally, the CAD system tries to mimic classical manufacturing approaches (such as turning, 
drilling, or milling) starting out from a form similar to a semimanufactured product. This makes 
the handling of the CAD tool easily accessible to engineers. However, this also makes it very 
complicated to use for isogeometric analysis. 

Sevilla, Fernandez-Mendes, and Huerta have proposed an alternative on middle ground between 
isogeometric analysis and standard finite elements: the NURBS-enhanced finite element method 
(NEFEM) 1119411196111951 1202 II . In NEEEM the boundary of the geometry is represented using 
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NURBS, thus leading to an exact representation, whereas for the interior, a standard finite element 
mesh, with all benehts of existing meshing algorithms, is utilized. Note that this leads to two 
different kinds of elements: (1) standard hnite elements in the interior and (2) elements with 
one NURBS edge alongside the boundary. Usually, elements of category (1) will be in the vast 
majority, keeping the computation very efficient. Through the elements of category (2), the 
exact geometry is made available during the process of evaluating the integrals of the appropriate 
weak form (e.g., of the governing equations Equations Q-Q): The integration domain O!^ is 
no longer only an approximation of the real domain Q.. The availability of the exact geometry is 
however not resorted to for the interpolation of the unknown solution. The latter is approximated 
using standard Lagrangian shape functions both along the boundary and in the interior. As a 
consequence, NEFEM is no longer an isoparametric method. This is in contrast to IGA. 



(a) Position of quadrature 
points 



(b) Example of a linear 
shape function 


Figure 27: (a) The NEFEM quadrature points are adapted to the curved triangle shape, (b) 
Example of a linear shape function in the NEFEM context. Negative values and values 
larger than 1 may occur. 

The main advantage of NEFEM over standard finite elements is that the position of the 
integration points needed in the finite element method are determined from the curved NURBS 
geometry and not from the approximated geometry (cf. Figure [T7|). If we consider boundary 
integral, such as the integral arising from the surface tension, the integration points are evenly 
distributed on the curved geometry. Due to the continued approximation of the unknown function 
through Lagrange polynomials an oddity arises: negative shape function values may occur along 
the boundary (cf. Figure [27]l. Furthermore, the shape function connected to all nodes, also those 
not located on the NURBS edge, are non-zero along the NURBS edge. This has the consequence, 
that interior nodes contribute to boundary integrals - a highly unusual situation in the finite 
element method. 

This section concludes the introduction of the numerical methods. The next sections will 
concentrate on standard applications. 


9 Bubbles and drops 

The simulation of drops (liquid) or bubbles (gas) in a surrounding bulk phase is one of the 
classical benchmark cases in multi-phase flow. It addresses the issues of surface tension effects, 
moving interfaces, as well as the question of break-up and coalescence of phase domains (in 
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this case, the drops). The test case is motivated hy the application of liquid-liquid or liquid-gas 
extraction columns Il28l . The relevant effects can however also he extended for example to the 
simulation of the alveoli in the lung. 


9.1 Static drop inside a rectangular domain 


The first test case regularly used in the context of drop simulations is a circular drop inside a 
rectangular domain (cf. Figure [28(a)| ). The domain Q .2 is occupied hy fluid 2, which is surrounded 
hy another fluid 1 residing in the domain fli. The dimensions of the domains differ within the 
literature. What remains in common is that the steady Stokes equations without external forces 
(i.e., in particular without gravity) are solved on the two domains (cf. Equation Q-Q). Density 
and viscosity of the two fluids are usually chosen to he equal. However, the presence of surface 
tension with a surface-tension coefficient of 7 is assumed. Since the analytical value of the 
sum of principal curvatures is known for circles {k = ^) and spheres (fc = ^), according to the 


Laplace-Young equations, the analytical solution is (cf. Section 2.1 ): 


u(x) =0, 


p{x) = 


0, X E , 

YK, X E D 2 . 


(91) 

(92) 



pix,y) 



(a) Computational (b) Exact pressure solution 

domain 


Figure 28: Static test case with circular interface. In (h), a zero reference pressure is assumed 
for fli. Furthermore, the parameters have been chosen to pi=p 2 = l.Okg/m^, pi=p 2 = 
l.Okgjmjs, 7 = l.Okg/s^ and r = 0.5m. 

With the analytical solution at hand, this test case is particularly valuable to demonstrate the 
capabilities of a chosen method to handle the jump of quantities across the interface as well as 
the evaluation of geometrical quantities of the interface. This is even more the case as Ganesan 
et al. Il 86 l point out that this solution is in general not represented well by a discretization 
scheme: the discrete velocities are not equal to 0 - we see spurious velocities. lf 86 l identifies two 
causes: inadequate representation of the pressure solution (particularly the pressure jump) and 
an insufficient approximation of the curvature. They give the error bound on as (in our own 
notation): 
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(93) 


|m |i <C 


( i 


inf 

62 * 


p-q llo + 


sup 


\{ k ’\ w ^ -n) — { k , w ^ -n) 


From theoretical considerations, Ganesan et al. 118611 conclude that the use of continuous 
approximation functions in the case of a discontinuous solution is to always he avoided - this 
holds true even if the discontinuities are resolved hy the mesh. The jump is then spread over several 
elements, thus distorting the solution (cf. Figure |29(a) I. Flere, and also in Gross and Reusken 
|[97l . it is shown that in general, the error in the pressure approximation inf^Agg/, \\ p — q’^ ||o is 
hounded hy cy/h, leading to very slow convergence. 



(a) Continuous interpolation. (b) Doubled interface nodes. 


Figure 29: Comparison of the pressure jump using continuous interpolation functions and inter¬ 
face nodes with doubled pressure degrees of freedom ||69]| . ©Wiley-Blackwell 

These theoretical results brought about a whole set of ideas aiming at numerically obtaining 
both an instant pressure jump and an improved curvature approximation. 

Pressure jump-. In the area of level-set methods, the work concentrated on developing problem- 
suitable finite element pressure spaces: XFEM evolved as the main means of incorporating 
the pressure jump (cf. Section [6^ . In ll97ll95l[T86]| . a Heaviside enrichment is employed for 
this purpose. In the case of continuous basis functions, convergence orders for the pressure of 
only 0{h^l^) in the L^-norm were reached. For the velocity, the reported convergence order is 
in the //'-norm. However, if the Heaviside-enrichment was employed in combination 
with an appropriate curvature approximation, a convergence order of with a > 1 could be 

recovered. Also in the area of level-set approximations, Ausas et al. lITlIl present an approach 
where the pressure space is locally enriched at those elements, which contain the interface. The 
corresponding elements are subdivided along the interface, and each part of the domain is then 
only influenced by pressure degrees of freedom on one side of the interface. The method does 
not introduce new degrees of freedom for the pressure. In 112011 . the stability of this method is 
compared to the method of Coppola-0 wen iHfil already mentioned in Section 

Along similar lines, but in the context of the PFEM method, the pressure degrees of freedom at 
the interface are duplicated, once associated to and once to Q. 2 , in 11581 . This leads to the 
possibility of an exact representation of the pressure jump even on coarse meshes. In the context 
of boundary-conforming meshes, the same idea has been employed in 1691 . The velocity values 
of the doubled nodes are coupled in the context of the boundary condition in Equation ( fTO] ), but 
the pressure is allowed to take on different values. A sample result is depicted in Eigurej^ The 
advantage of the local approaches is that they take features of the individual problem into account, 
e.g., the knowledge that the pressure jump will only ever occur at the interface. This makes the 
enrichment less general, but usually much more efficient and less detrimental to the numerical 
scheme (such as, for example, the ill-conditioned matrices XFEM can induce). 
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Curvature approximation: Most approaches never compute the curvature of the interface, but 
resort to the Laplace-Beltrami technique or the continuum surface force technique as described in 
Section (among many others lITSl '). Among the rare explicit approaches, several researchers 
have attempted the use of splines as an interface representation, from which the curvature can then 
be derived analytically. This includes ll^ where piecewise interpolated cubic splines provide an 


) represents the entire circle, thus providing an excellent means for evaluating the curvature at 
any given point of the interface. 12021 furthermore employs a space-time version of the NEFEM 
(cf. Section [O] ). Here, a combination of PIP 1 finite elements and a quadratic NURBS for the 
geometry description is employed. With this combination, the analytical solution can be obtained 
exactly even on a coarse mesh with only 8 EE nodes on the interface. 


additional description of the discretized interface. In 11691120211 . a single NURBS (cf. Section 8.4 


9.2 Rising drops or bubbles 

The test case described in the previous section can be extended when independent material 
parameters (/ii ,2 and pi 2 ) as well as gravity are taken into account. Depending on the density 
relation, the drop in D 2 will then either rise or fall in due to buoyancy. In comparison to the 
previous test case, the set-up of this section needs to additionally consider the velocity kink at the 
interface, the jump in the pressure gradient, and most of all, a scheme for deforming domains. 

Experimentally, the phenomenon of a rising drop is for example described in ll45ll . The 
book defines three dimensionless numbers, any two of which completely characterize the flow 
regime. They involve the gravity g, the drop diameter d, viscosity p and density p, as well as the 
surface-tension coefficient 7 : 


• The Reynolds number, relating inertia to friction forces. 


Re = 


P2y/gdd 




(94) 


• the Edtvds number, relating buoyancy forces and surface tension. 


Eo = 



(95) 


• and the Morton number, relating viscous forces and the surface tension. 


Mo = 


P27^ 


Eo'^ 

R^' 


(96) 


Depending on these characteristic numbers, three regimes can be distinguished: spherical, 
ellipsoidal, and spherical cap. Spherical drops occur at low Re and Eo. The ellipsoidal regime is 
located around relatively high Re and intermediate Eo, while both high Re and Eo yield a drop in 
the spherical cap regime. 

Hysing et al. published benchmark computations for a single rising bubble 111 1611 . It contains 
two benchmark cases, one in the ellipsoidal regime and one in the spherical cap regime, featuring 
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1.0 



Figure 30: Rising bubble: Initial computational domain. 


thin filaments and break-up. The results of three different numerical methods are described; two 
use the finite element method in conjunction with level-set and one finite elements in an ALE 
setting. 

In rtl9ll . the first benchmark of HI 161 was repeated. Here, the surface tension effects dominate 
the bubble shape and no break-up should occur. The rectangular domain has the size of 1.0 x 2.0m 
with an initially circular bubble with diameter d = 0.5m, as shown in Figure [30| 

The properties of the fluids are p\ = lOOkg/m^, p 2 = 1000%/m^, pi = \kg/m/s, pi = 
X^kgjmjs, fy = —g = —0.9Sm/s^, and surface tension coefficient 7 = 24.5kg/s^. Under these 
circumstances, the characteristic Reynolds and Ebtvbs numbers can then be evaluated to: 


Ml 


(97) 


Eo = ^ = 10. 


(98) 


The spatial resolution is 80 x 160 elements and At = 0.002^. No-slip boundary conditions are 
assumed at the top and bottom boundary, slip boundary conditions are used along the vertical 
walls, and zero pressure is specified at the upper boundary. As an initial condition, tbe velocity 
field is set to 0 . Tbe results - Figure 31(a) reports tbe bubble shape over time - give further 
support to the conclusions in II116II . In Figure 31(a)[ the bubble shape at t = 3.0^ is compared 
with the result from II116I . showing a very good agreement. In order to compare a time-evolving 
quantity, the rise velocity of the bubble, defined by the authors of the benchmark as: 


_IaiV{x,t)dx 
^rise — r 1 j 

In, ^dx 

is depicted in Figure [3l(b)| 


( 99 ) 
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(a) Bubble shape (b) Rise velocity 


XFEM OTP2D 


Figure 31: Rising bubble: Comparison of the bubble shape and rise velocity obtained with the 
XFEM level-set approach in II191II and simulation data from II116II . 


Other publications where this or comparable benchmarks have been studied are, e.g., the 
following. In 111581 . a particle finite element method (PFEM) is employed to evaluate the second 
test case. In this method, the marker particles are treated in a Lagrangian manner. The interaction 
between the particles is portrayed through a connecting mesh, which is renewed after each 
displacement update. In addition, the article discusses the difficulties in modeling break-up and 
coalescence, which are always mesh-depended. Aland et al. consider both test cases with a phase 
field method ||4|. The governing equations for the mixture of incompressible, immiscible fluids 
lead to a coupling of the momentum equation ([T]) with a Cahn-Hillard type phase field equation. 
These are solved using an adapfive finife elemenf scheme. Also wifh fhe phase-field mefhod, rising 
bubble examples similar fo II116II are invesfigafed in ll98l . A fhermodynamically consisfenf sef of 
governing equafions is presented. The numerical scheme combines finife elemenf discrefizafions 
for the momentum, phase-field, and chemical potential equations with a finite volume approach 
for all convective terms. Very large deformations as well as hreak-up and coalescence are inherent 
for these methods. In II129I the authors were able to reproduce both benchmark cases (with and 
without breakup) with the finite volume method. Qualitatively, the results show the same features 
as in mH; the quantified results, however, present significant deviations from the ones obtained 
with the finite element codes. Doyeux et al. utilize the benchmark to validate their finite element 
method with interface capturing using level-set before they move on to the simulation of vesicles 
(a structure which serves as an imitation of the mechanical behavior of red blood cells). Again, 
both benchmark cases are presented ll62ll . In ||3l, the benchmark case is extended to 3D. The 
results of three different codes: XFEM with level-set, finite differences with level-set, and finite 
volumes with VOE are compared. 
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10 Sloshing tank 


A second common application in free-surface flows are sloshing tanks. The numerical simulation 
of sloshing tank problems is motivated by, e.g., the analysis of storage tanks under seismic loading, 
or partially filled tanks on ships such as fuel tanks and maritime transport of fluids or liquid 
gas. The aim of the simulation is to gain information about the forces acting on the surrounding 
structures, e.g., the tank walls or the ship, in order to investigate possible failure mechanisms. 
When considering such effects, the phenomenon of sloshing is in general not negligible when 
investigating the responses of the tank structure to the load due to the large mass of the liquid 
in combination with the low stiffness of the tank walls II174I . If the problem is to be considered 
to its full extent, it falls into the category of fluid-structure interaction problems, where the thin 
structure of the tank wall reacts to the forces generated by the motion of the fluid and vice-versa. 
In this work, we will concentrate on the fluid component of such simulations. 

The amplitude of the sloshing depends on the amplitude and frequency of the seismic load as 
well as the fluid fill level, the fluid properties, and the tank geometry II174I . In the next section, we 
will focus on the last point, including arbitrary tank geometries. Usually, the considered tanks are 
either rectangular or cylindrical. With regard to the fluid-structure-interaction problems, methods 
for arbitrarily shaped tanks need to be developed, however, in order to incorporate the different 
failure mechanisms of the tank wall into the sloshing phenomenon. Known shapes under failure 
of the tank structure are elephant footing, diamond-shaped buckling, as well as random convex 
and concave bumps II1281I . 

Within the overall context of the fluid-structure-interaction problem, concessions are often 
made with regard to the accuracy of the fluid flow model. Many authors decide to model the liquid 
in the tank as inviscid, incompressible and irrotational. This enables the use of a simple Laplace 
equation, which describes the velocity potential. If the full Navier-Stokes equations ([T]l-Q are 
chosen as governing equations, the sloshing tank is an application where the body forces are 
of great importance. In general, a gravitational force will be considered, since otherwise, any 
upward sloshing could never be reversed. In cases of rotating tanks, Coriolis forces are added 
to the general body forces. If seismic loading is to be considered, an alternating acceleration 
(usually orthogonal to gravity) is imposed. 

Another factor, which categorizes different scenarios, is the regularity of the flow conditions. 
While moderate deformations can be accurately and efficiently handled by interface tracking 
methods, wave breaking, air entrapment, and touching of the tank ceiling is more easily conducted 
with interface-capturing methods. Benchmark problems in this area are in general restricted 
to rectangular tanks with either an initially flat surface, which is subsequently excited (e.g.. 
Ill 10112081 ). or a free surface with a prescribed original deformation, which is then allowed to 
settle into equilibrium (e.g., in iTMllSlI l. Even though this is a two-phase scenario, with a liquid 
phase in the container secluded by a gaseous phase on the top, in this application, the modelling 
is usually restricted to the liquid phase, as the density of the gaseous phase is several orders of 
magnitude lower ||92l . Exceptions need to be made, e.g., if the flow conditions are so violent 
that gas entrapments occurs. Ansari et al. lEl use the modal method, a reduced-order method, to 
model tank sloshing based on the flow potential (Laplace equation). The authors investigate the 
question under which conditions the top fluid can be neglected during the simulation. The results 
were later further analyzed in ||9^ . In ll88l . the scope of the modal method is extended from 
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vertical walls to tapered walls. Numerical simulation and experimental validation of sloshing 
in a 2D rectangular tank is presented in |[53l . The full Navier-Stokes equations are solved and 
the interface is captured using marker particles, which are updated in a Lagrangian manner. 
Subsequently, a global mass correction scheme is applied in order to ensure conservation of 
mass. Sloshing liquid motion in combination with a floating roof, which underlies the non-linear 
equation of motion, is the topic of II153L The flow is modelled through the analytical solution 
modes of the velocity potential, and only the interface to the floating roof requires finite-element 
discretizations. In II174L the authors decide to solve the compressible Navier-Stokes equations 
with an ALE finite element method. 

In a similar application, Lee et al. II133II simulate non-linear free-surface flows around ship 
bows including wave breaking phenomena. Due to the wave breaking, interface-capturing 
methods are more appropriate then interface tracking. In the referenced paper, a modified marker- 
density method is employed, circumventing the oversimplifications a traditional MAC-method 
would involve. As a solution method, finite differences with a sequential solution of velocity and 
pressure are used. One important challenge in this type of simulation is the conjunction point 
between the ship hull and the free surface, namely imposing the appropriate boundary conditions. 

10.1 The boundary conditions 

The free surface requires a no-penetration boundary condition (cf. Equation (fTO])), indicating that 
the shape of the free surface is dictated through the fluid velocity. If one would like to consider 
surface tension effects in the pressure solution. Equation ( [TT| ) needs to be incorporated. Note 
that any technique relying on partial integration now needs to take the boundary integral into 
account (in contrast to the situation described in Section]^, as the free surface in a sloshing tank 
is not closed. An exception can be made if we assume that the contact angle 6, i.e., the angle 
between the tangent vector along E'"', n', and the boundary F is assumed to be constant at 90° (cf. 
Figurej^. Then the boundary term is again 0 as shown in |[84l . This generalization is acceptable 
if wetting effects, as described in II125L are assumed to not play an influential role. This holds 
true for the sloshing tank: As Behr points out in l[22l . the capillary effects are negligible in the 
case of a sloshing tank, thus precluding the need for a specific wetting model. Instead, a global 
slip condition is found appropriate to allow the free surface to slide freely along the tank walls. 
The slip boundary condition can be expressed in the following way ll2^ : 

u-n = 0 onr^//p 

t[C7(u,p)n = 0 on F.up (100) 

t2<r(u,p)n = 0 onr.iip. 

As an alternative, a Navier slip condition II197I could be employed to incorporate at least some 
amount of wall friction. 

The tangential and normal vectors can either be computed analytically, if an appropriate 
description is available, or from the finite element mesh. 
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Figure 32: Close-up of the contact line F^ between an interface F'”' and a domain boundary F. 

The tangent vector of the interface at F^, n', and the boundary F span the contact angle 
a. If wetting models are not considered, a = 90°. 

10.2 Tanks with walls of arbitrary shape 

If interface-capturing methods are used, the shape of the tank wall is insignificant, as long as the 
computed velocity field complies with the slip boundary condition. Since the marker particles, 
the level-set field, or the VOF field are advected with the fluid velocity, they will obey the shape 
of the boundary. In an interface-tracking context, tanks with walls of arbitrary shape have so far 
posed a significant challenge. The reasons for this are two-fold. One reason lies in the mesh 
deformation, which needs to comply with the analytic boundary at all times. This means that in 
addition to the discretized boundary, sufficient information on the analytic boundary needs to be 
stored, e.g., in form of a CAD model. The second reason is due to the discrepancy between the 
discrete slip condition and the analytic slip condition. In Il22l . it was found that this remains true 
even if a conforming normal vector is utilized. 

The general procedure for incorporating a slip boundary condition into an interface-tracking 
context is to first assemble the full system matrix A, irrespective of any boundary conditions 
that need to be applied. Subsequently, all equations are rotated in such a way, that they are now 
formulated in a coordinate system spanned by the tangent and the normal vector of the respective 
boundary. Only then are the Dirichlet conditions, e.g., zero velocity in normal direction, applied. 
This again requires two steps: (1) At any finite element node (or degree of freedom) with a 
Dirichlet boundary condition, the corresponding weighting function is equal to zero. Therefore, 
the corresponding equation is deleted from the matrix A. (2) Furthermore, all occurrences of the 
specified degree of freedom are moved to the right-hand-side. Note that, usually, we encounter 
homogeneous Dirichlet conditions, where nothing needs to be done in this second step. The 
described procedure is straightforward if the tangential coordinate system is aligned with the 
main coordinate axes. In all other cases, the tangential coordinate system together with the 
corresponding rotation matrix remains to be determined. 

Consider the original linear system of equations with the nodal velocities in the (x,y,z)- 
coordinate system: 
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( 101 ) 


(All Ai2 
A21 A22 

V : 



For each node i with degrees of freedom u; = Ui^y a rotation matrix O, is 

determined, which transforms from the (x, 3 ^, 2 )-coordinate system to a given tangential coordinate 
system. Rotation matrices are square, orthogonal matrices with determinant 1 ( otherwise 
O^* = would not hold and the rotation could not he easily reversed). They can he obtained 
hy placing the maps of the basis vectors of W‘^‘^ into the columns of the rotation matrix. The 
local rotation matrix then takes the form HTl : 


0,= (ti t 2 n). (102) 

With the help of the rotation matrix, new degrees of freedom u, are defined in the tangential 
coordinate system: 


U; = OiUi 


(103) 


Equation (101 1 may then be formulated in terms of the new degrees of freedom: 


With 


/All Ai2 

A21 A22 

V ^ ^ 


/Oi 0 

0 O 2 

V: 


we can solve the system using: 



/Oi 

0 ...\ 


/or' 0 ...\ 


0 

0 

0 

V ^ 

02 


1 tN ... 

0 

0 • • • 


... 0 

... 0 


(104) 


(105) 


/ 0 [ 0 


0 o[ 


V 


/Oi 0 

0 O 2 

V: : 



/o[ 0 
c 

V: 


0 of 



(106) 


For analytical domain shapes (straight or tilted lines, circles), 11221 already shows how the 
rotation matrix O can be formed. In order to cope with arbitrary domain shapes, we propose a 
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boundary description that is based on NURBS (cf. Section 8.4 for the introduction to NURBS 
and the notation). The use of NURBS as a boundary description has the following advantages: 
(1) NURBS can represent arbitrary analytical and free-form shapes, (2) they can be extracted 
from CAD-models, (3) tangent and normal vectors can be analytically computed at any point on 
the NURBS. Due to point (3), the alignment of the system of equations of the flow solution can 
be performed in the standard way using a rotation matrix O obtained as in Equation ( 102 ). 

The first tangent vector is obtained as the first derivative of the NURBS surface (or curve) 
with respect to the first local parameter d ; a second tangent vector can be obtained by taking the 
derivative with respect to the second local parameter i: 


U,NURBS 


dsie,i) 
de ’ 


h,NURBS 


ds{e,i) 

di 


(107) 


Note, however, that in physical space, these two tangent vectors are not orthogonal on each 
other. This only holds for the parameter space spanned by 6 and i. Since the rotation matrix O 
requires orthogonality in physical space, Gram-Schmidt orthogonalization is performed to obtain 
the two tangent vectors i\ „ and i 2 .n- 


^ _ h NURBS 

— II . II » 

II NURBS II 

, _ ^2NURBS — {^2NURBS-^l,n)h,n 


^2NURBS — {^2NURBS ' 


(108) 

(109) 


From the two tangent vectors, the normal vector can be computed via cross product. 

Applying a similar procedure to the mesh deformation scheme is slightly more complex. The 
aim is that the boundary nodes of the mesh slide along the arbitrarily shaped boundary in order to 
accommodate the deformation imposed through the free surface. The individual displacement for 
each node 1) is to be determined via the mesh deformation equation ( |7^ . However, as for the 
flow solution, we require a slip boundary condition. In addition to this condition, which restrict 
the nodal displacement to the local tangent direction, the requirement that the nodes remain on 
the NURBS at all times needs to be fulfilled. In total, the following steps are performed: 


1. Rotate the mesh equation into the local tangential coordinate systems. 

2. Drop the equation responsible for the normal component. 

3. Compute the tangential displacement. 

4. From the displacement in x-coordinates, deduce a displacement on the NURBS in the local 
NURBS coordinates 6 and i. 


The steps 1. through 3. are performed within the solution of the mesh deformation equation. 
Step 4. involves a transformation of the displacement given in tangential coordinates into the 
local NURBS coordinates. If the displacement is directly available in NURBS coordinates, it is 
ensured that all nodes remain on the NURBS. The following relation holds for the displacement 
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from one iteration (it) to the next (it + 1), where an iteration can be both a linearization step (such 
as in the Newton-Raphson algorithm) or a time step: 




-x,v=S(0,,+i)-S(0,v) ^ JA0. 

linearization 


( 110 ) 


Ad refers to the displacement in NURBS-coordinates and J is the Jacobian of S with respect 
to 6. From Equation ( |110| ), A6 can be computed. 

Figure [^illustrates the resulting displacement. 



Figure 33: The nodal displacement l) is computed in terms of the tangential coordinate system 
spanned by ty„ , , n. This displacement is transformed into a displacement expressed 

in the NURBS parametric coordinates 6 based on a local linearization of the NURBS 
curve. 


In 2D, J = If and Ad can be readily obtained. In 3D, cannot be determined easily and 
the (overdetermined) equation system needs to be solved. Starting point is step 3., where the 
tangential displacement Ati , At 2 has been obtained. The displacement reads 


(tl,n 



( 111 ) 


This displacement can be exactly represented as a displacement in the ti,t 2 -coordinate system 
- i.e., the coordinate system aligned with the NURBS local coordinates. What remains is to 
determine the appropriate magnitude of the displacement, A6,Ai, which fulfils the relation 



( 112 ) 


Note that (ti 


12 ) is also the Jacobian matrix J. The above relation can be rewritten as 


At\t\^n + At2t2,n = AFt] + Alt2 . 


(113) 


With the aim of comparing coefficients, the expressions (1091 are inserted in Equation (113l 
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A0tj “t" Ait2. 


(114) 


All 


■+Al 2 - 


h — (t2 •tl^n)tl,« 
h — (t 2 


leading to the final result of 


AG 

At 


tl II V II t2 — (t2 11/ 

A?2 

t2 — (t2 II 


(115) 

(116) 


In such cases, one final remark needs to he made about the intersection between the free surface 
and the wall (F^ in Figure [32]) . Here, it is important that the displacement of the free surface is 
chosen compatible with the displacement along the wall. 

The effect of the new implementation of the slip condition is detailed by means of a numerical 
example. 

In the test case, we consider an undeformable tank with curved side walls as indicated in 
Figure]^ The tank wall is defined through a quadratic NURBS curve with 18 control points. 
Along all walls, a slip boundary condition is imposed. The top boundary constitutes a free-surface. 
The tank is subjected to a downward gravitational acceleration g = —0,9Sm/s^ and a horizontal 
sinusoidal acceleration s = —0.2sin(t). The considered fluid has a density of p = lOOOkg/m^ 
and a viscosity of p = OAkgjmjs. Figure depicts the velocity vectors of the flow solution 
after 4^ and 13.5^. Note how the velocity aligns with the boundary. The mesh deformation is 
conform with the boundary at all times. 


11 Conclusion 

This paper gave an overview of the challenges in the simulation of fluid flow in deformable do¬ 
mains. Five main methods have been introduced with their individual advantages and drawbacks. 
As was seen, some may be used more often than others, but each has its niche. A recent area for 
the Marker and Cell method has been found in flow simulations for movies, where a qualitatively 
correct solution is more important than quantitative accuracy. Level-set has its main advantage 
when topological changes are to be expected within the domain. The volume-of-fluid method 
is utilized frequently in commercial codes. In cases where boundary conditions are difficult to 
impose, the use of the phase-field method can be worthwhile. The interface tracking approach 
may have a limited scope of applications, but comes with very high computational efficiency. 

Out of a very large number of possibilities, the applications of drops and sloshing tanks were 
selected to demonstrate some of the approaches when computing free-boundary problems. 
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(b) t = 13.5J 

Figure 34: Sloshing in a tank with curved side walls: the tank wall is described hy a 
quadratic NURBS curve with knot vector [0 0 0 1/16 2/16 3/16 4/16 5/16 6/16 
7/16 8/16 9/16 10/16 11/16 12/16 13/16 14/16 15/16 1 1 1] and control points 
(10,10);(10,5);(10,5);(12,4);(8,3);(12,2);(8,1);(10,0);(10,0);(0,0);(0,0);(2,1); 
(—2,2); (2,3); (—2,4); (0,5); (0,5); (0,10). The velocity vectors of the flow solution 
are depicted after 4^ and 13.55. Note how the velocity aligns with the boundary. 12001 
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